{ "cells": [ { "cell_type": "markdown", "id": "117fc616-fd5c-45c3-a910-3003f656bf0b", "metadata": {}, "source": [ "## 6 - Precision in Predictions\n", "### Figuring out precision and other types of tests" ] }, { "cell_type": "markdown", "id": "3efc371a-93a3-4a11-a27a-4607cf222138", "metadata": {}, "source": [ "There are only a few examples this week to drive home the idea of precision in predictions, or equivalently, how *un*certain we are of our predictions.\n", "\n", "The main take home message is that our predictions will vary in precision, but that the impact of this precision varies with the magnitude of the effect. A large prediction that is imprecise may be of more use than a tiny, precise prediction - as ever, it all depends on the use of the model. Traditionally, p-values will be significant when precision is high, with comes with larger sample sizes. We must not be misled by significance, and should think carefully about what kind of hypotheses are important about our predictions." ] }, { "cell_type": "markdown", "id": "75cde913-c26c-4959-a6fc-821de92d48c3", "metadata": {}, "source": [ "### Finding precision\n", "We will see some examples of precision in prediction. First, lets import all the usual packages we need:" ] }, { "cell_type": "code", "execution_count": 1, "id": "d62fb6a0-1b28-4f65-8c8e-290096251043", "metadata": {}, "outputs": [], "source": [ "# Import what we need\n", "import pandas as pd # dataframes\n", "import seaborn as sns # plots\n", "import statsmodels.formula.api as smf # Models\n", "import marginaleffects as me # marginal effects\n", "import numpy as np # numpy for some functions\n", "\n", "# Set the style for plots\n", "sns.set_style('whitegrid')\n", "sns.set_context('talk')" ] }, { "cell_type": "markdown", "id": "293e0a10-0fdd-4667-9b8b-957791e5161b", "metadata": {}, "source": [ "## Intervention data\n", "We will load the data for the childhood intervention study to examine prediction and precision in some detail, and explore how to set specific hypotheses to test.\n", "\n", "We download it from here: https://vincentarelbundock.github.io/Rdatasets/csv/mlmRev/Early.csv" ] }, { "cell_type": "code", "execution_count": 2, "id": "1737dd45-3ea2-4e03-97de-0d502f88c4c9", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
rownamesidcogagetrt
01681031.0Y
12681191.5Y
2368962.0Y
34701061.0Y
45701071.5Y
\n", "
" ], "text/plain": [ " rownames id cog age trt\n", "0 1 68 103 1.0 Y\n", "1 2 68 119 1.5 Y\n", "2 3 68 96 2.0 Y\n", "3 4 70 106 1.0 Y\n", "4 5 70 107 1.5 Y" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Read in \n", "cog = pd.read_csv('https://vincentarelbundock.github.io/Rdatasets/csv/mlmRev/Early.csv')\n", "cog.head()" ] }, { "cell_type": "markdown", "id": "b04ffc0c-3b1a-4aa1-ab47-d4432fbe22bc", "metadata": {}, "source": [ "And fit a simple model to the data:" ] }, { "cell_type": "code", "execution_count": 3, "id": "49075383-c40d-43f2-b4f4-14ad2c767e44", "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
OLS Regression Results
Dep. Variable: cog R-squared: 0.314
Model: OLS Adj. R-squared: 0.309
No. Observations: 309 F-statistic: 69.97
Covariance Type: nonrobust Prob (F-statistic): 9.45e-26
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
coef std err t P>|t| [0.025 0.975]
Intercept 124.5216 2.950 42.206 0.000 118.716 130.327
trt[T.Y] 9.4903 1.497 6.340 0.000 6.545 12.436
age -18.1650 1.819 -9.988 0.000 -21.744 -14.586


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified." ], "text/latex": [ "\\begin{center}\n", "\\begin{tabular}{lclc}\n", "\\toprule\n", "\\textbf{Dep. Variable:} & cog & \\textbf{ R-squared: } & 0.314 \\\\\n", "\\textbf{Model:} & OLS & \\textbf{ Adj. R-squared: } & 0.309 \\\\\n", "\\textbf{No. Observations:} & 309 & \\textbf{ F-statistic: } & 69.97 \\\\\n", "\\textbf{Covariance Type:} & nonrobust & \\textbf{ Prob (F-statistic):} & 9.45e-26 \\\\\n", "\\bottomrule\n", "\\end{tabular}\n", "\\begin{tabular}{lcccccc}\n", " & \\textbf{coef} & \\textbf{std err} & \\textbf{t} & \\textbf{P$> |$t$|$} & \\textbf{[0.025} & \\textbf{0.975]} \\\\\n", "\\midrule\n", "\\textbf{Intercept} & 124.5216 & 2.950 & 42.206 & 0.000 & 118.716 & 130.327 \\\\\n", "\\textbf{trt[T.Y]} & 9.4903 & 1.497 & 6.340 & 0.000 & 6.545 & 12.436 \\\\\n", "\\textbf{age} & -18.1650 & 1.819 & -9.988 & 0.000 & -21.744 & -14.586 \\\\\n", "\\bottomrule\n", "\\end{tabular}\n", "%\\caption{OLS Regression Results}\n", "\\end{center}\n", "\n", "Notes: \\newline\n", " [1] Standard Errors assume that the covariance matrix of the errors is correctly specified." ], "text/plain": [ "\n", "\"\"\"\n", " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: cog R-squared: 0.314\n", "Model: OLS Adj. R-squared: 0.309\n", "No. Observations: 309 F-statistic: 69.97\n", "Covariance Type: nonrobust Prob (F-statistic): 9.45e-26\n", "==============================================================================\n", " coef std err t P>|t| [0.025 0.975]\n", "------------------------------------------------------------------------------\n", "Intercept 124.5216 2.950 42.206 0.000 118.716 130.327\n", "trt[T.Y] 9.4903 1.497 6.340 0.000 6.545 12.436\n", "age -18.1650 1.819 -9.988 0.000 -21.744 -14.586\n", "==============================================================================\n", "\n", "Notes:\n", "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n", "\"\"\"" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Model\n", "cog_mod = smf.ols('cog ~ age + trt', data=cog).fit()\n", "cog_mod.summary(slim=True)" ] }, { "cell_type": "markdown", "id": "1d893787-150a-43ff-b6fe-ec7e603c97e2", "metadata": {}, "source": [ "We can then predict the effect of the intervention for 2 year olds easily:" ] }, { "cell_type": "code", "execution_count": 4, "id": "9c0655c8-b7bd-4514-9e89-d8f03c40087e", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "shape: (2, 13)
trtagerowidestimatestd_errorstatisticp_values_valueconf_lowconf_highrownamesidcog
stri64i32f64f64f64f64f64f64f64i64i64i64
"Y"2097.6818441.34388472.686230.0inf95.04788100.3158073112096
"N"2188.191551.44528961.0199920.0inf85.35883591.0242653112096
" ], "text/plain": [ "shape: (2, 9)\n", "┌─────┬─────┬──────────┬───────────┬───┬─────────┬─────┬──────┬───────┐\n", "│ trt ┆ age ┆ Estimate ┆ Std.Error ┆ … ┆ P(>|z|) ┆ S ┆ 2.5% ┆ 97.5% │\n", "│ --- ┆ --- ┆ --- ┆ --- ┆ ┆ --- ┆ --- ┆ --- ┆ --- │\n", "│ str ┆ str ┆ str ┆ str ┆ ┆ str ┆ str ┆ str ┆ str │\n", "╞═════╪═════╪══════════╪═══════════╪═══╪═════════╪═════╪══════╪═══════╡\n", "│ Y ┆ 2 ┆ 97.7 ┆ 1.34 ┆ … ┆ 0 ┆ inf ┆ 95 ┆ 100 │\n", "│ N ┆ 2 ┆ 88.2 ┆ 1.45 ┆ … ┆ 0 ┆ inf ┆ 85.4 ┆ 91 │\n", "└─────┴─────┴──────────┴───────────┴───┴─────────┴─────┴──────┴───────┘\n", "\n", "Columns: trt, age, rowid, estimate, std_error, statistic, p_value, s_value, conf_low, conf_high, rownames, id, cog" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Predict for 2 year olds\n", "p = me.predictions(cog_mod,\n", " newdata=me.datagrid(cog_mod,\n", " trt=['Y', 'N'],\n", " age=2)\n", " )\n", "\n", "p" ] }, { "cell_type": "markdown", "id": "452cafce-49cb-441e-bd28-c8ebc698a04c", "metadata": {}, "source": [ "The `conf_low` and `conf_high` represent the 95% confidence interval bounds. " ] }, { "cell_type": "code", "execution_count": 5, "id": "8c727974-64be-4dd7-80fb-a4b32ec8a8d3", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmsAAAHHCAYAAADgeh/sAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAABHhUlEQVR4nO3deXRURcL+8aeTdJZOgBgwEMISE1mHIKAYWRReZWAcRVnEkU1kF0TFV4YRRVGBAU5QXCCK4IAwQGYUXhdQEQIu7JtAgCwSloRlCERCgtnT/fuDX9qJSUiHbsg1/f2c45nTt+pW1dU5J8+pW1XXZLPZbAIAAIAheVT3AAAAAFAxwhoAAICBEdYAAAAMjLAGAABgYIQ1AAAAAyOsAQAAGBhhDQAAwMC8qnsAcE58fLwKCwvl4eEhHx+f6h4OAABwQH5+vqxWq8xmsyIjI69al7D2O1dYWCibzabi4mLl5ORU93AAAEAVFBYWVlqHsPY75+HhoeLiYplMJvn5+VX3cAAAgANyc3Nls9nk4VH5ijTC2u+cj4+PcnJy5Ofnp1atWlX3cAAAgAMSEhKUk5Pj0BImNhgAAAAYGGENAADAwAhrAAAABkZYAwAAMDDCGgAAgIER1gAAAAyMsAYAAGBghDUAAAADI6wBAAAYGGENAADAwAhrAAAABkZYAwAAMDDCGgAAgIER1gAAAAzMq7oHcC1OnDihPn36qF+/fnrllVfKrbNt2zYtWrRIiYmJysvLU3h4uB577DE98sgjMplMZeoXFRVp9erV+te//qUTJ07Iy8tLkZGRGjNmjKKioqo8xi+//FLLli3TsWPHVFxcrJYtW2rYsGHq2bNnlduqDnmFxXrv2xT773HdI+Rr9qzGEQEA4J5+d2HtwoULGj9+vHJzcyuss2rVKr366qsym82KioqS2WzWjh07NHXqVO3bt0+zZs0qVd9ms+lvf/ub1q5dqzp16qhz587KzMzUtm3btHXrVs2cOVP9+/d3eIzR0dFavHixLBaLoqKiVFBQoF27dunpp5/WhAkT9PTTT1/z898o+UVWvR33k/33iK63ENYAAKgGv6uwlpCQoGeffVYnT56ssM7x48c1ffp0BQQEaPny5WrdurUk6cyZMxo2bJjWrFmjbt266U9/+pP9ntWrV2vt2rVq1aqVli5dqsDAQEnS9u3bNXbsWL322mvq3LmzQkJCKh3j9u3btXjxYoWEhGjFihUKDQ2VJCUmJuqJJ57QggUL1K1bN7Vt29aJfxMAAMBd/C7WrF26dEnR0dF69NFHdfLkSTVq1KjCuosXL1ZxcbFGjhxpD2qS1LBhQ/sr08WLF5e6Z+HChZKkqVOn2oOaJHXq1EnDhg1Tfn6+li9f7tBYS9p67rnn7EFNklq2bKmJEyfKZrPpww8/dKgtAACA30VYW7ZsmRYvXqygoCC999576tOnT4V1N2/eLEnlrg3r3LmzatWqpfj4eJ07d06SdPToUaWmpqpu3bq6/fbby9zTq1cvSVJcXFyl47x8+bJ27dolT09P3XfffWXKe/bsKZPJpG+//VbFxcWVtgcAAPC7eA3aoEED/e1vf9OgQYPk6+urw4cPl1vvwoULysjIkNlsVnh4eJlyT09PhYeH68CBA0pKSlL9+vWVnJwsSWrWrFm5Gw9Krqempio3N1d+fn4VjjMlJUXFxcVq3LixAgICypQHBQWpbt26unDhgk6cOKGIiAhH/xVUymazuTQAWn/TlrW4WMXFv4tsDwCA4dlsNofr/i7C2oABAxyql56eLkmqV6+ePDzKDxbBwcGl6pb8b/369cut7+Pjo9q1a+vSpUu6cOGCGjduXGn/FbVV0v+FCxeUnp7u0rCWm5ur/fv3u6y9XwqspX7Hx8fL35uwBgDAjVaj/vrm5ORIknx9fSus4+PjU6ruL7/84vA9JXUrUpW2SvoHAAC4mt/FzJqjSmbTynud+Vsl04+enp4O31OZqrRVlelPR/j5+alFixYuay8rt1D67Nd1epGRkartZ3ZZ+wAAuLOkpKSrHkP232pUWPP395ck5eXlVVgnPz9fkmSxWKp8z9XWq11r/65iMpnsYdEVPDytv/nt6dL2AQBwZ1WZJKpRr0FL1opduHChwpmrknVlJWvXSu45f/58ufXz8/OVlZUlk8mkm2++2aH+K2qrvP4BAACupkaFtcDAQNWvX18FBQXlHpxbXFysY8eOSZL9lWHJ//70009l6pdct9lsaty4caWzYbfeeqs8PT11+vTpcqc2f/75Z2VkZMjHx0dNmzat0rMBAAD3VKPCmiR169ZNkvTNN9+UKdu6dauys7PVsmVLNWjQQJLUtGlThYWFKT09vdzdlF9//bUkqXv37pX27ePjo7vuukuFhYXatGlTmfL169fLZrOpa9euMptZ/wUAACpX48La4MGD5enpqYULF5YKX2fOnNH06dMlSWPHji11z+OPPy5Jevnll3XhwgX79e3bt2vZsmUym80aMWJEqXvS09OVkpJif63527bmzJmjEydO2K8nJibq7bffLrd/AACAitSoDQbSr591euONNzRo0CDdeeed8vX11c6dO5WTk6MBAwboz3/+c6l7Bg4cqO+//17ffvutevXqpaioKGVnZ2vPnj2y2WyaPXt2me+Cvvnmm/q///s/9e3bV7Nnz7Zf7969ux577DHFxsbqoYce0l133aXi4mLt3LlThYWFevbZZ3XbbbfdkH8XAADg96/GhTVJGjNmjCIiIrR06VIdPHhQJpNJERERGjhwoPr27VumvoeHh+bPn69//vOfWrNmjbZs2aKAgAB17txZY8eO1Z133lml/l999VW1bdtWq1at0q5du+Tj46PbbrtNw4cPV48ePVz1mAAAwA2YbK4+8As3VEJCgnJycmSxWNSqVSuXtXspt1C3vfbrur8D03qqDuesAQDgElX5+13j1qwBAADUJIQ1AAAAAyOsAQAAGBhhDQAAwMAIawAAAAZGWAMAADAwwhoAAICBEdYAAAAMjLAGAABgYIQ1AAAAAyOsAQAAGBhhDQAAwMAIawAAAAZGWAMAADAwwhoAAICBEdYAAAAMjLAGAABgYIQ1AAAAAyOsAQAAGBhhDQAAwMAIawAAAAZGWAMAADAwwhoAAICBEdYAAAAMjLAGAABgYIQ1AAAAAyOsAQAAGBhhDQAAwMAIawAAAAZGWAMAADAwwhoAAICBEdYAAAAMjLAGAABgYIQ1lJGelafl20+UurZ8+wmlZ+VVz4AAAHBjXtU9ABhH/KlLev/7FK0/9B8VWW2lyuZ+k6y3Nv6kXm0a6Ml7IhTZqE41jRIAAPdCWIMk6Z87TuqVzw7pNxmtlCKrTesOntVX8Wf1+sNtNOSupjdugAAAuCnCGvTPHSc19dNDDte32qSpnx6SySQNjiKwAQBwPbFmzc3Fn7qkVz5zPKj9t5c/PaT4U5dcPCIAAPDfCGtu7v3vU6766vNqrDZp4fcprh0QAAAopca+Bv3qq6+0cuVKHTlyRDabTWFhYerfv78GDBggb29ve71Tp07pvvvuc6jNCRMm6Omnn6603pYtWzRy5MgKyy0Wi3788UeH+rye0rPytP7Qf5xq4+tD/1F6Vp6Ca/u6aFQAAOC/1ciwNm3aNMXGxkqSmjdvrkaNGikxMVGvv/66vvzyS8XExKhOnSu7GS0Wi3r37l1hW5mZmfrhhx8kSa1atXKo/yNHjkiSIiMjFRYWVqbcx8enKo9z3Xx+4EyZXZ9VVWS16fMDZzTq7nAXjQoAAPy3GhfWPvvsM8XGxspsNis6Olr333+/JKmwsFBz5szR8uXLNWPGDEVHR0uSgoKCNHfu3HLbslqtGjVqlCTpySefVI8ePRwaw6FDV9aATZw4UV27dnX2ka6blPOXXdTOLy5pBwAAlFXj1qytWrVKkjR69Gh7UJMks9msF154QREREfr888+VnJxcaVsxMTHaunWr2rVrp2eeecbhMZTMrLVp06aKo7+xfskvdlE7RS5pBwAAlFXjwlpSUpIklTsL5uXlpY4dO0qSvvvuu6u2c+LECb3//vsym82aOXOmPD09Heo/KytLaWlpatKkiQIDA6s2+BvM38exZ6q8nRo3QQsAgGHUuL+yxcVXZosCAgLKLffyuvLIKSlX38U4ffp0FRYWasSIEbr11lsd7v/w4cOSpCZNmmjBggX6+uuvlZqaqoCAAHXq1Enjx49XeLgx1ndF3Fz+v6Oqt+PvknYAAEBZNS6shYeHKyEhQbt27VLTpqUPbLXZbNq3b58kKSMjo8I2du3apS1btshisWj06NFV6r8krG3ZskW7d+9Wx44dFRISosOHD+uLL75QXFycYmJi1KlTpyo+2dXZbDZ7UHXUA23qa/ZXiU5tMvDyMOmBNvWr3DcAAO7MZnP8b2+NC2v9+vXTzJkzNXfuXDVr1kzt2rWTdGWzwLvvvmtfT1ZQUFBhG4sXL5Yk/eUvf1FQUFCV+i9pv2PHjnrrrbdUr149e3+zZ8/WihUrNHHiRG3YsEG1a9eu6uNVKDc3V/v376/yfXc29NG2U9f+gfaoUB+dOZaoM9fcAgAAuBqTrSrR7neguLhYTz/9tOLi4uTh4aHIyEjVq1dPSUlJOnfunPr376/Y2Fh17dpVH374YZn7U1JS9MADD8jLy0ubNm1ScHBwlfovKCjQqVOnFBwcXOZVbHFxsfr376+EhARNnTpVQ4cOdepZJSkhIUE5OTnXfH/KxUK9sDFD1mu410PS7B51FXGT+Zr7BwDAnVkslkqPBqtxM2uenp6aP3++Vq5cqY8//lgJCQmyWCyKiorS/Pnzdfz4ccXGxtrPWfuttWvXymazqUuXLlUOapLk7e1d4Zo0T09Pde/eXQkJCYqPj69y21fj5+enFi1aVPm+dpLyLal6+fMjVb73tYdbq/+dTap8HwAA7i4pKUm5ubkO1a1xYU2SPDw8NGTIEA0ZMqRM2YYNGyRJoaGh5d67fv16SdKDDz54XcYWEhIiSQ7/B3KUyWRyeMfqbw3tfIs8PD308qeHHPr0lIdJmt6nDR9xBwDgGplMJofr1rijO1JTU7VlyxadP3++3PJt27ZJktq2bVumLC0tTSkpKTKbzQ5/guq/5efn66WXXtK4ceN08eLFcuucPXtWktSgQYMqt389DY5qqs+e6qoH24bIy6P8/wN5eZj0YNsQffZUV4IaAAA3SI0La6tXr9bIkSP1ySeflCk7cuSI9u/fr8DAQHXp0qVM+YEDByRJf/jDH2SxWKrct4+Pj7Zu3apNmzYpLi6uTHlBQYHWrVsnSerevXuV27/eIhvV0fxBHbTthXs1qWfzUmWTejbXtin3av6gDopsVP4rZAAA4Ho1Lqz16NFDJpNJS5cuVVpamv362bNnNWnSJNlsNo0dO7bcMFayjqxkB+nVZGdnKyUlRampqaWuDxo0SJIUHR2txMRE+/W8vDxNmTJFqamp6tixY7lh0SiCa/tqaKewUteGdgpTcC0+1g4AwI1W49asRUZGatSoUVq0aJF69+5t/2LBzp07lZ+fr759++qJJ54o996ScNekSeWL5jds2KApU6YoNDRUmzZtsl8fPny49u3bp82bN6t///7q0KGDAgMDtXfvXmVkZCg8PFzz5s1z/kEBAIBbqHFhTZKef/55NW7cWKtWrdKOHTvk7++vtm3batCgQbr//vsrXNT3888/S3JuPZnZbFZMTIw+/vhjrV69WocOHZLValXjxo01aNAgjRgx4ppesQIAAPdU485Zczcl56w5ck5LVVzKLdRtr31j/31gWk/V8eM8NQAAXKEqf79r3Jo1AACAmoSwBgAAYGCENQAAAAMjrAEAABgYYQ0AAMDACGsAAAAGRlgDAAAwMMIaAACAgRHWAAAADIywBgAAYGCENQAAAAMjrAEAABgYYQ0AAMDACGsAAAAGRlgDAAAwMMIaAACAgRHWAAAADIywBgAAYGCENQAAAAMjrAEAABgYYQ0AAMDACGsAAAAGRlgDAAAwMMIaAACAgRHWAAAADIywBgAAYGCENQAAAAMjrAEAABgYYQ0AAMDACGsAAAAGRlgDAAAwMMIaAACAgRHWAAAADIywBgAAYGCENQAAAAMjrAEAABgYYQ0AAMDACGsAAAAGRlgDAAAwMK/qHsD18tVXX2nlypU6cuSIbDabwsLC1L9/fw0YMEDe3t5l6v/rX//SK6+8UmF7zZo109q1ax3uPz4+XjExMTp06JCysrLUuHFj9enTR8OGDZPZbL6mZwIAAO6nRoa1adOmKTY2VpLUvHlzNWrUSImJiXr99df15ZdfKiYmRnXq1Cl1z+HDhyVJUVFRCg4OLtNmSEiIw/1v3rxZEyZMkNVq1R133KHatWtr9+7dio6O1rZt27Rw4UICGwAAcEiNC2ufffaZYmNjZTabFR0drfvvv1+SVFhYqDlz5mj58uWaMWOGoqOjS91XEtZee+013XLLLdfc/6VLlzRp0iRJ0qJFi9S1a1dJUmZmpsaMGaOtW7dq2bJlGjly5DX3AQAA3EeNW7O2atUqSdLo0aPtQU2SzGazXnjhBUVEROjzzz9XcnKyvaywsFDJycmqVauWwsLCnOp/xYoVunz5svr06WMPapIUGBioWbNmSZKWLFmi4uJip/oBAADuocaFtaSkJElSjx49ypR5eXmpY8eOkqTvvvvOfv3o0aMqKChQmzZtZDKZnOp/06ZNkqSePXuWKYuIiFDz5s11/vx5HTx40Kl+AACAe6hxYa1kxiogIKDcci+vK29+U1JS7NdKXoHWr19fc+bM0Z/+9Ce1bdtW3bp107Rp03Tu3DmH+//pp58kXVkrV55bb71VkpSYmOhwmwAAwH3VuDVr4eHhSkhI0K5du9S0adNSZTabTfv27ZMkZWRk2K+XhLVPP/1UAQEBuuOOOxQSEqLDhw8rNjZWGzZs0JIlS9SiRYur9n3p0iXl5eVJuhL8ylOyeSE9Pf3aHrACNpvNpa9Wrb9py1pcrOLiGpftAQCoFjabzeG6NS6s9evXTzNnztTcuXPVrFkztWvXTpJktVr17rvv6siRI5KkgoIC+z0l13r16qW///3v9lm57OxsvfTSS1q/fr2eeeYZrVu3zj4zV56cnBxJkre3tzw8yg82vr6+peq6Sm5urvbv3++y9gqKbXq0tb/9d+KRQ/L2dO4VMQAAqLoaF9YGDx6sHTt2KC4uTgMHDlRkZKTq1aunpKQknTt3To899phiY2NLha6PPvpIaWlpatq0aakz2GrVqqVZs2bpxx9/1IkTJ/T999/r3nvvrbDvkoDmyLq3qiTq6uDtadJf/lCruocBAIDbq3FhzdPTU/Pnz9fKlSv18ccfKyEhQRaLRVFRUZo/f76OHz+u2NjYUues+fr6qlmzZuW25+/vr7vuukuff/654uPjrxrW/P2vzETl5+fLarWWO7tW8prUz8/Pmccsw8/Pr9LXtAAAwBiSkpKUm5vrUN0aF9akKzNcQ4YM0ZAhQ8qUbdiwQZIUGhrqcHslB+JW9i81ICBAAQEBunz5ss6fP1/uurWStWrlHbzrDJPJJE9PT5e2CQAAro+qnD5R41aMp6amasuWLTp//ny55du2bZMktW3bVpJ07tw5TZkyRc8880yFC/TPnj0rSWrQoEGl/ZfsAj169Gi55SXXmQUDAACOcFlYu3Tpkr7++mvFxMRozpw5kq68Dty1a5erunDI6tWrNXLkSH3yySdlyo4cOaL9+/crMDBQXbp0kXRlXdq6deu0fv36cseamZmpzZs3y2Qy6Z577qm0/27dukmS1q9fX6YsJSVFycnJCgoK0m233VbVRwMAAG7I6bBms9n0zjvvqHv37nruuef07rvvaunSpZKkU6dOadiwYRo4cKB+/vlnZ7tySI8ePWQymbR06VKlpaXZr589e1aTJk2SzWbT2LFjZbFYJEkWi0V9+vSRJL366qs6ffq0/Z5Lly7p2WefVXZ2th566CGFh4fby7Kzs5WSkqLU1NRS/ffr108BAQFavXq14uLi7NczMzP14osvSpJGjBjBt0EBAIBDTDYntyX+9a9/1dq1a2Wz2XTTTTcpPz9fubm5SkhI0N69ezV48GCZTCaFh4frk08+cfnC+vLMnTtXixYtkp+fn/2LBTt37lR+fr769u2rv//976UW/1++fFnDhw/XwYMH5evrqw4dOsjX11e7d+9Wdna2br/9di1atMi+gUCS1qxZoylTpig0NNT+1YISX3zxhSZPniybzab27durbt262r17tzIzM3X33Xfrvffec1lYS0hIUE5OjiwWi1q1auWSNgEAwPVVlb/fTs2sffPNN/riiy8UFBSkRYsWafv27WrZsqW9/Pbbb9eKFStUt25dHTt2TMuWLXOmO4c9//zzev311xUWFqYdO3YoPj5ebdu21bx58zRr1qwyuzQDAgK0YsUKTZo0SWFhYdq7d6927NihJk2aaMqUKfroo49KBbXK9O7dW8uXL9fdd9+to0ePauvWrapfv76mTJmimJgYZtUAAIDDnJpZGzlypLZt26YlS5borrvukiQNGjRIP/74oxISEuz19uzZoyFDhqh169Zas2aN86OGHTNrAAD8/tywmbVDhw4pJCTEHtQqcscddyg0NFQnTpxwpjsAAAC341RYy8nJUWBgoEN1g4KCVFRU5Ex3AAAAbsepsFavXj2dPHmy0k8nFRYW6sSJE6pXr54z3QEAALgdp8LanXfeqZycHK1cufKq9T766CNlZ2frjjvucKY7AAAAt+NUWBs+fLg8PDw0Z84cLVu2TBcvXixVnpGRobfeektvvvmm/RNQAAAAcJzT56ytWLFCM2bMKHM9MDBQmZmZkq4cnPvss89q3LhxznSFcrAbFACA358bthtUkgYPHqyFCxeqRYsWstls9n8uXrwom82mJk2a6M033ySoAQAAXAMvVzRyzz336J577tHp06f1008/KTs7W35+frrlllsUERHhii4AAADcklNhbffu3apVq5b9qwWhoaEKDQ0tt+6WLVt08uRJDR482JkuAQAA3IpTr0GHDh1a7nq18rz11lt66623nOkOAADA7Tg8s3b58uUyuz0lKS8vT2lpaRXeZ7PZdPr0aR07dqzS89gAAABQmsNh7ZdfflHv3r2Vn59vv2YymXT48GH17NnToTbatWtX5QECAAC4M4dfg9avX1/Dhw8vteNTUqnfV/snJCREU6dOvW4PAgAAUBNVaYPB+PHj9cgjj0i6EtJ69OihyMjIq65F8/DwkMViUZ06dZwaKAAAgDuqUlgzm82ldnt27NhRLVq0qHAHKAAAAJzj1NEdy5cvd9U4AAAAUA6nv2DwW1artdQ/BQUFysrKUnJyshYuXOjq7gAAAGo0p79g8P333+vdd99VcnKyCgoKKq0/duxYZ7sEAABwG06FtUOHDmn8+PEqLi6u9Aw1Ly8vtW/f3pnuAAAA3I5TYW3ZsmUqKipSs2bNNGrUKPn6+urZZ59Vr1699Je//EX/+c9/tGbNGu3Zs0d33HGHli5d6qJhAwAAuAenwtqePXvk6empd955R7fccoskKSQkRGlpaercubMkqW/fvnrmmWe0ceNGrVu3Tg888IDzowYAAHATTm0wyMjIUMOGDe1BTZJatmxZav2ayWTSlClTJElr1qxxpjsAAAC34/Ru0MDAwFK/w8LCVFxcrOPHj9uvNWzYUE2bNlVSUpKz3QEAALgVp8Ja3bp1df78+VLXGjVqJEk6evRoqev+/v7KzMx0pjsAAAC341RYa9Omjc6dO6cffvjBfi08PFw2m027d++2X8vLy1NqaiqfnAIAAKgip8Laww8/LJvNpmeeeUbR0dEqKipS+/btVadOHX3yySf69NNPlZycrKlTpyo7O7vU2jYAAABUzqmwdt999+mBBx5Qbm6uPvroI3l6esrX11cDBw5UUVGRpkyZoocffljr1q2TyWTS448/7qpxAwAAuAWnv2DwxhtvqEuXLtq6datMJpMkacKECTp37pw+/fRT2Ww2eXp6avjw4erZs6fTAwYAAHAnJltlnx5wQnp6us6cOaMmTZooKCjoenXj1hISEpSTkyOLxaJWrVpV93AAAIADqvL32+mZtasJDg5WcHDw9ewCAACgRnNJWCsoKNCxY8d0+fLlSr8R2rFjR1d0CQAA4BacDmvz5s3TsmXLlJeXV2ldk8mkI0eOONslAACA23AqrC1dulQLFy50uP51XB4HAABQIzkV1v7973/LZDKpX79+euqppxQcHCwvr+u6DA4AAMCtOJWs0tLSVLduXU2fPl0eHk5/ZhQAAAC/4VTCCggI0M0330xQAwAAuE6cSlkdO3bU8ePH9csvv7hqPAAAAPgvToW18ePHq7i4WK+99hqbBwAAAK4Dp9astWzZUm+//bYmTJig/fv3q0uXLgoKCrJ/dqo8EyZMcKZLAAAAt+JUWLt8+bI++OADWa1WpaWlKTY2ttJ7blRY++qrr7Ry5UodOXJENptNYWFh6t+/vwYMGCBvb+8y9XNzc7VkyRKtX79eJ0+elNVqVWhoqO677z6NGjVKgYGBDve9ZcsWjRw5ssJyi8WiH3/88VoeCwAAuBmnwtq8efO0f/9+SdLNN9+shg0bymw2u2JcTpk2bZo9ODZv3lyNGjVSYmKiXn/9dX355ZeKiYlRnTp17PUzMzM1dOhQJScnq1atWmrfvr08PT0VHx+vRYsW6csvv9TKlSvVoEEDh/ovOfg3MjJSYWFhZcp9fHycf0gAAOAWnAprcXFxMplM+vvf/66+ffu6akxO+eyzzxQbGyuz2azo6Gjdf//9kqTCwkLNmTNHy5cv14wZMxQdHW2/Jzo6WsnJybrzzjv19ttv2z86n5WVpf/93//VDz/8oJdfflmLFi1yaAyHDh2SJE2cOFFdu3Z18RMCAAB34tQGg59//llNmzY1TFCTpFWrVkmSRo8ebQ9qkmQ2m/XCCy8oIiJCn3/+uZKTkyVJeXl5Wrt2rSRp9uzZ9qAmSbVr19acOXNkMpn0ww8/6OLFiw6NoWRmrU2bNi55JgAA4L6cCmvBwcHy9PR01VhcIikpSZLUo0ePMmVeXl72D8l/9913kqSMjAz94Q9/UIcOHRQaGlrmnrp166pOnTqy2WxKT0+vtP+srCylpaWpSZMmVVrnBgAAUB6nXoP26tVL//jHP3T48GH94Q9/cNWYnFJcXCzpyoG95Sn5HFZKSookKTQ0VCtXrqywvZMnTyozM1MeHh6qX79+pf0fPnxYktSkSRMtWLBAX3/9tVJTUxUQEKBOnTpp/PjxCg8Pr9IzAQAA9+VUWBs3bpw2btyocePG6cUXX9Tdd98tf39/V43tmoSHhyshIUG7du1S06ZNS5XZbDbt27dP0pUZNUfMnTtXktS5c2eHZspKwtqWLVu0e/dudezYUSEhITp8+LC++OILxcXFKSYmRp06darCU1XOZrPZgyoAADC2qpxP61RYe/3119W4cWNt2bJFzz33nEwmk2rVqiU/P79y65tMJm3evNmZLivVr18/zZw5U3PnzlWzZs3Url07SZLVatW7775rX09WUFBQaVvvv/++vvnmG/n6+mry5MkO9V/SfseOHfXWW2+pXr169v5mz56tFStWaOLEidqwYYNq1659DU9YvtzcXPvOXAAAUHOYbE58eqBly5ZV68xkUkJCwrV255Di4mI9/fTTiouLk4eHhyIjI1WvXj0lJSXp3Llz6t+/v2JjY9W1a1d9+OGHFbbzzjvvaMGCBfLw8NDcuXP1wAMPONR/QUGBTp06peDg4DKvYouLi9W/f38lJCRo6tSpGjp0qFPPKkkJCQnKyclxuh0AAHDjWSwWtWrV6qp1nJpZmzVrljO3Xxeenp6aP3++Vq5cqY8//lgJCQmyWCyKiorS/Pnzdfz4ccXGxpY6Z+2/FRQUaOrUqfrss8/k5eWlWbNmORzUJMnb27vCNWmenp7q3r27EhISFB8ff03PVxE/Pz+1aNHCpW0CAIDrIykpSbm5uQ7VdSqsGenIjv/m4eGhIUOGaMiQIWXKNmzYIEnl7vzMyMjQU089pR9//FEBAQF6++23XX5OWkhIiCQ5/B/IUSaTyXA7cwEAQPmu9mnO33Lq6A4jSk1N1ZYtW3T+/Plyy7dt2yZJatu2bZn7BgwYoB9//FGNGjWyvyqtivz8fL300ksaN25chWeynT17VpIc/hoCAABwbw7PrG3fvl2S1KFDB/vnkkquVYWrd0H+1urVq/X+++9r4sSJGjduXKmyI0eOaP/+/QoMDFSXLl3s18+dO6fHH39cZ8+eVdu2bfX++++rbt26Ve7bx8dHW7du1dmzZxUXF6dHHnmkVHlBQYHWrVsnSerevXvVHw4AALgdh8Pa8OHD5eHhoXXr1umWW26xX6vKNJ7JZLLvlrxeevTooYULF2rp0qV68MEH1bhxY0lXZrQmTZokm82msWPHymKx2O/561//qrNnz6pFixZaunSpQ8ePZGdnKz09XWazWU2aNLFfHzRokN544w1FR0erTZs29k0YeXl5eumll5SamqqOHTuWCosAAAAVqdKaNavVWuZaVTaTOrHx1GGRkZEaNWqUFi1apN69e9u/WLBz507l5+erb9++euKJJ+z1t27dqp07d0qSatWqpWnTplXY9rPPPmsPfxs2bNCUKVMUGhqqTZs22esMHz5c+/bt0+bNm9W/f3916NBBgYGB2rt3rzIyMhQeHq558+ZdhycHAAA1kcNhLTEx0aFrRvD888+rcePGWrVqlXbs2CF/f3+1bdtWgwYN0v33319qNvC/z33bs2fPVdsdNmyYPaxVxGw2KyYmRh9//LFWr16tQ4cOyWq1qnHjxho0aJBGjBhRalYPAADgapw6Z62qMjIyrmktGCpWcs6aI+e0AAAAY6jK32+ndoPed999eu655xyq+9hjjxn2qA8AAACjciqsnT59Wunp6ZXWs1qtOn/+fIXHWQAAAKB8Dq9ZO3r0aLmL75OTkzV48OAK77PZbDp37pzOnDmjhg0bXtsoAQAA3JTDYe3WW2+Vr6+vtm7dar9mMpmUnZ2tvXv3OtRGeV8UAAAAQMWqdHTHyy+/rLVr19p/z58/Xw0bNlS/fv0qvMdkMsnf31+tWrVSVFTUtY8UAADADVUprIWFhWnChAn23/Pnz1dISEipawAAAHAdpz7kHhcXZ//0FAAAAFzPqbAWGhpaYVleXp62bdsmq9WqO+64Q4GBgc50BQAA4JacCmvSlY+gv/fee2rYsKHGjBkjSUpJSdHw4cN1/vx5SZKfn59mzJihP//5z852BwAA4FacCms///yzHn30UaWnp6t79+7266+88orS09PtmwsuX76syZMnq0WLFoqIiHB2zAAAAG7DqUNxP/roI507d05NmjTRX/7yF0nSyZMntXfvXnl6emrVqlXas2ePxowZo6KiIi1dutQVYwYAAHAbToW177//Xl5eXvrwww/tM2vffvutJKlDhw5q166dJOnpp59W7dq1tWPHDme6AwAAcDtOhbW0tDSFhYWpUaNG9mvbtm2TyWRS586d7dfMZrMaNWrk0KepAAAA8CunwlpxcbG8vb3tv4uKirR7925J0p133lmqbm5urkwmkzPdAQAAuB2nwlpoaKhOnz6twsJCSdLu3buVk5Mjf39/+ytQ6cqO0bS0NIWEhDg1WAAAAHfjVFiLjIxUVlaW5s6dq8TERL311lsymUzq1q2bPD09JUkZGRn661//quLiYnXq1MklgwYAAHAXToW10aNHy9fXV8uWLVPfvn114MABeXp6avTo0ZKkPXv2qFu3btq9e7dq1aqlESNGuGTQAAAA7sKpsBYeHq5//OMfioyMlLe3t5o3b6733ntPLVu2lCQFBwerqKhIzZo106pVq0ptRAAAAEDlnP6CQfv27fXvf/+73LJGjRrp008/tYc3AAAAVI3TYa3EDz/8oM2bN+vYsWPKzs7W6tWrdfnyZW3cuFH169fXTTfd5KquAAAA3IbTYS0jI0MTJ07Unj17JEk2m81+RMeZM2c0f/58LV++XB988IFuu+02Z7sDAABwK06tWSsoKNDIkSO1e/du+fv7649//KPq16//a+MeHgoMDNSlS5c0fPhwnT592ukBAwAAuBOnwtqKFSuUmJiodu3a6ZtvvtE777yj0NBQe3nz5s21ceNGtW/fXrm5uVqyZInTAwYAAHAnToW1devWycPDQ9HR0QoKCiq3TkBAgObOnStPT0/98MMPznQHAADgdpxas3bs2DFFRESocePGV60XGhqqsLAwpaamOtMdAADADZNXWKz3vk2x/x7XPUK+Zs8bPg6nwprVanW4rtlstn/VAAAAwOjyi6x6O+4n++8RXW+plrDm9LdBT5w4ocuXL1+13sWLF/XTTz+VWs8GAACAyjkV1rp166bCwkJFR0dftd6MGTNUXFysu+++25nuAAAA3I5Tr0FHjhyp1atX69///rcyMjLUu3dvZWdnS5JSUlKUnJysFStWaO/evfL399cTTzzhijEDAAC4DafCWt26dRUTE6Px48dr48aNiouLs5c9+OCDkq4ckmuxWPTmm2+WOoMNAAAAlXPqNagk3X777fr888/1+OOPKyQkRDabzf5P3bp19cgjj+jTTz/VPffc44rxAgAAuBWXfBu0fv36evHFF/Xiiy8qJydH2dnZslgsqlWrliuaBwAAcFsu+5B7CYvFIovF4upmAQAA3JLTr0EBAABw/RDWAAAADIywBgAAYGCENQAAAAMjrAEAABgYYQ0AAMDACGsAAAAG5vJz1oziq6++0sqVK3XkyBHZbDaFhYWpf//+GjBggLy9vcu9Z9u2bVq0aJESExOVl5en8PBwPfbYY3rkkUdkMpmq1H98fLxiYmJ06NAhZWVlqXHjxurTp4+GDRsms9nsikcEAABuoEbOrE2bNk0TJ07Url271LBhQ0VFRenixYt6/fXXNXz4cF26dKnMPatWrdLw4cO1e/dutW7dWlFRUUpJSdHUqVP14osvVqn/zZs367HHHtO3336rsLAwde3aVenp6YqOjtbYsWNVWFjoqkcFAAA1XI2bWfvss88UGxsrs9ms6Oho3X///ZKkwsJCzZkzR8uXL9eMGTMUHR1tv+f48eOaPn26AgICtHz5crVu3VqSdObMGQ0bNkxr1qxRt27d9Kc//anS/i9duqRJkyZJkhYtWqSuXbtKkjIzMzVmzBht3bpVy5Yt08iRI1396AAAoAaqcTNrq1atkiSNHj3aHtQkyWw264UXXlBERIQ+//xzJScn28sWL16s4uJijRw50h7UJKlhw4Z65ZVX7HUcsWLFCl2+fFl9+vSxBzVJCgwM1KxZsyRJS5YsUXFx8bU/JAAAcBs1LqwlJSVJknr06FGmzMvLSx07dpQkfffdd/brmzdvliT17NmzzD2dO3dWrVq1FB8fr3PnzlXa/6ZNmypsKyIiQs2bN9f58+d18OBBB54GAAC4uxr3GrRkxiogIKDcci+vK4+ckpIiSbpw4YIyMjJkNpsVHh5epr6np6fCw8N14MABJSUlqX79+lft/6effpIkNW/evNzyW2+9VcnJyUpMTFT79u0deygH2Gw2ZusAAHAh62/+rlqLi1Vc7Jp5LpvN5nDdGhfWwsPDlZCQoF27dqlp06alymw2m/bt2ydJysjIkCSlp6dLkurVqycPj/L/AwQHB5eqW5FLly4pLy9PkioMdY62VVW5ubnav3+/S9sEAMCd/VJgLfU7Pj5e/t43/qVkjXsN2q9fP0nS3LlzS4UXq9Wqd955R0eOHJEkFRQUSJJycnIkSb6+vhW26ePjU6puRUrKvb29Kwx+Jf1U1hYAAIBUA2fWBg8erB07diguLk4DBw5UZGSk6tWrp6SkJJ07d06PPfaYYmNj7a9DS0KVI+eoVTZl6cq2qsrPz08tWrRwaZsAALizrNxC6bM4++/IyEjV9nPNWalJSUnKzc11qG6NC2uenp6aP3++Vq5cqY8//lgJCQmyWCyKiorS/Pnzdfz4ccXGxqpOnTqSJH9/f0myv74sT35+viTJYrFcte+StvLz82W1WsudXSvpx8/Pr+oPdxUmk0menp4ubRMAAHfm4Wn9zW9Pl/2trcph+zUurElXZriGDBmiIUOGlCnbsGGDJCk0NFTSr2vLLly4IJvNVu6/vJL1ZSXrzSoSEBCggIAAXb58WefPny933ZqjbQEAAEg1cM1aamqqtmzZovPnz5dbvm3bNklS27ZtJV05/6x+/foqKCjQyZMny9QvLi7WsWPHJMmh14wlu0CPHj1abnnJdV5ZAgAAR9S4sLZ69WqNHDlSn3zySZmyI0eOaP/+/QoMDFSXLl3s17t16yZJ+uabb8rcs3XrVmVnZ6tly5Zq0KBBpf2XtLV+/foyZSkpKUpOTlZQUJBuu+02h58JAAC4rxoX1nr06CGTyaSlS5cqLS3Nfv3s2bOaNGmSbDabxo4dW2r92eDBg+Xp6amFCxeW2kF65swZTZ8+XZI0duzYUv1kZ2crJSVFqamppa7369dPAQEBWr16teLifl2UmJmZaf/G6IgRI/iYOwAAcEiNW7MWGRmpUaNGadGiRerdu7f9iwU7d+5Ufn6++vbtqyeeeKLUPS1bttTEiRP1xhtvaNCgQbrzzjvl6+urnTt3KicnRwMGDNCf//znUvds2LBBU6ZMUWhoqP2rBdKVtWivvvqqJk+erKeeekrt27dX3bp1tXv3bmVmZuruu+8u0z8AAEBFalxYk6Tnn39ejRs31qpVq7Rjxw75+/urbdu2GjRokO6///5yNxGMGTNGERERWrp0qQ4ePCiTyaSIiAgNHDhQffv2rVL/vXv3VkhIiH2mLjExUY0bN9a4ceM0aNAgZtUAAIDDTDZXH/iFGyohIUE5OTmyWCxq1apVdQ8HAIAa41JuoW577df17Aem9VQdF52zVpW/3zVuzRoAAEBNQlgDAAAwMMIaAACAgRHWAAAADIywBgAAYGCENQAAAAMjrAEAABgYYQ0AAMDACGsAAAAGRlgDAAAwMMIaAACAgRHWAAAADIywBgAAYGCENQAAAAMjrAEAABgYYQ0AAMDACGsAAAAGRlgDAAAwMMIaAACAgRHWAAAADIywBgAAYGCENQAAAAMjrAEAAPxGelaelm8/Uera8u0nlJ6Vd8PH4nXDewQAADCo+FOX9P73KVp/6D8qstpKlc39JllvbfxJvdo00JP3RCiyUZ0bMibCGgAAgKR/7jipVz47pN9ktFKKrDatO3hWX8Wf1esPt9GQu5pe93ER1gAAgNv7546TmvrpIYfrW23S1E8PyWSSBkdd38DGmjUAAODW4k9d0iufOR7U/tvLnx5S/KlLLh5RaYQ1AADg1t7/PuWqrz6vxmqTFn6f4toB/QZhDQAAuK30rDytP/Qfp9r4+tB/rusuUcIaAABwW58fOFNm12dVFVlt+vzAGReNqCzCGgAAcFsp5y+7qJ1fXNJOeQhrAADAbf2SX+yidopc0k55CGsAAMBt+ft4uqid63caGmENAAC4rYibA1zUjr9L2ikPYQ0AALith25rKC8Pk1NteHmY9NBtDV00orIIawAAwG0F1/ZVrzYNnGrjT20aKLi2r4tGVBZhDQAAuLUn74nQtU6ueZiksfdEuHZAv+3jurYOAABgcJGN6uj1h9tc073T+7RRZKM6Lh5RaYQ1AADg9obc1VQz+7ZxeIbNwyTN7Nvmun/EXZKu3z5TA/juu++0dOlSxcfHKy8vT/Xr11f37t01btw41atXT5J06tQp3XfffQ61N2HCBD399NOV1tuyZYtGjhxZYbnFYtGPP/7o2EMAAIAbYnBUU7UNDdTC71P09aH/lPtlAy8Pk/7UpoHG3hNx3WfU7H3ekF6qwZIlSzR79myZTCa1b99eQUFBOnjwoP75z39q/fr1WrlypZo0aSKLxaLevXtX2E5mZqZ++OEHSVKrVq0c6vvIkSOSpMjISIWFhZUp9/HxqfoDAQCA6y6yUR3NH9RB6Vl5+veeNM39JtleNqlncz3asbGCa12/zQTlqZFh7fTp03rjjTfk7e2txYsXKyoqSpJUUFCgyZMn66uvvtLMmTO1cOFCBQUFae7cueW2Y7VaNWrUKEnSk08+qR49ejjU/6FDhyRJEydOVNeuXV3wRAAA4EYKru2roZ3CSoW1oZ3CVMfPfMPHUiPXrG3fvl2FhYXq0qWLPahJkre3tyZOnChJ2rlzZ6XtxMTEaOvWrWrXrp2eeeYZh/svmVlr0+baFisCAACUqJFhzdPzyqcj0tPTy5RduHBBknTTTTddtY0TJ07o/fffl9ls1syZM+1tViYrK0tpaWlq0qSJAgMDqzZwAACA36iRr0E7deoks9msw4cPa9q0aRozZoyCgoK0f/9+vfLKK5KkMWPGXLWN6dOnq7CwUCNGjNCtt97qcN+HDx+WJDVp0kQLFizQ119/rdTUVAUEBKhTp04aP368wsPDr/3hKmCz2VRc7JqP0QIAAMn6m7+r1uJiFRe7Zp7LZiu7eaEiJltVav+ObNiwQS+99JIuXbpU6vpNN92k6dOn649//GOF9+7atUtDhw6VxWJRXFycgoKCHO538eLFio6OlnRlI0HHjh3l6empw4cP68KFC7JYLIqJiVGnTp2u7cF+IyEhQTk5OS5pCwAA/OqXAqse/+zXt3TLHg6Wv7drX0paLJZKNzDWyNegktS6dWv17NlTXl5eat++vf7nf/5HwcHBunjxoj744AOlpaVVeO/ixYslSX/5y1+qFNSkX9erdezYUZs2bdKHH36oDz74QJs3b9bgwYOVk5OjiRMnKisr69ofDgAAuI0aObOWkJCg4cOHy8fHR++9955at24tSSosLNSbb76pf/zjHwoJCdGXX34pi8VS6t6UlBQ98MAD8vLy0qZNmxQcHFylvgsKCnTq1CkFBwcrICCgVFlxcbH69++vhIQETZ06VUOHDnXuQfXrzJqfn59atGjhdHsAAOCKrNxCtZ8RZ//949T7VNtFu0GTkpKUm5vr0MxajVyzNmPGDF28eFELFiywBzVJMpvNmjx5sg4cOKC9e/fqk08+0eOPP17q3rVr18pms6lLly5VDmrSlR2nFa1J8/T0VPfu3ZWQkKD4+Pgqt301JpPJ4U0QAACgch6e1t/89nTZ31qTyfGPkda416D5+fnat2+fTCaTunTpUqbcZDKpW7dukn49D+2/rV+/XpL04IMPXpfxhYSESJJyc3OvS/sAAKBmqXFhLSsrS1ar9aozTSXXi4qKSl1PS0tTSkqKzGazw5+g+m/5+fl66aWXNG7cOF28eLHcOmfPnpUkNWjQoMrtAwAA91PjwlrdunUVGBgoq9Wqb7/9ttw6W7dulVT281EHDhyQJP3hD38os5bNET4+Ptq6das2bdqkuLi4MuUFBQVat26dJKl79+5Vbh8AALifGhfWPDw8NHDgQEnSzJkzlZz862cirFar5s+fr23btql27drq379/qXtL1pG1a9eu0n6ys7OVkpKi1NTUUtcHDRokSYqOjlZiYqL9el5enqZMmaLU1FR17Nix3Fe0AAAAv1UjNxg89dRTSkxM1ObNm/Xwww+rQ4cOqlOnjhITE3X69GlZLBa9/fbbZY7lKDnOo0mTJpX2sWHDBk2ZMkWhoaHatGmT/frw4cO1b98+bd68Wf3791eHDh0UGBiovXv3KiMjQ+Hh4Zo3b55rHxgAANRYNTKsmc1mvffee1qzZo3WrFmjxMRE5efnKzg4WI8++qhGjx5dbiD7+eefJTm3nsxsNismJkYff/yxVq9erUOHDslqtapx48YaNGiQRowYcU2vWAEAgHuqkeesuZOSc9YcOacFAAA47lJuoW577Rv77wPTeqqOi85Zq8rf7xq3Zg0AAKAmIawBAAAYGGENAADAwAhrAAAABkZYAwAAMDDCGgAAgIER1gAAAAyMsAYAAGBghDUAAAADI6wBAAAYGGENAADAwAhrAAAABkZYAwAAMDDCGgAAgIER1gAAAAyMsAYAAGBghDUAAAADI6wBAAAYGGENAADAwAhrAAAABkZYAwAAMDDCGgAAgIER1gAAAAyMsAYAAGBghDUAAAADI6wBAAAYGGENAADAwAhrAAAABkZYAwAAMDDCGgAAgIER1gAAAAyMsAYAAGBghDUAAAADI6wBAAAYGGENAADAwAhrAAAABkZYAwAAMDDCGgAAgIF5VfcAAAAAjMjHy0PP3tes1O/qUKPD2nfffaelS5cqPj5eeXl5ql+/vrp3765x48apXr16per+61//0iuvvFJhW82aNdPatWsd7js+Pl4xMTE6dOiQsrKy1LhxY/Xp00fDhg2T2Wy+5mcCAAA3hq/ZU8/9sXl1D6PmhrUlS5Zo9uzZMplMat++vYKCgnTw4EH985//1Pr167Vy5Uo1adLEXv/w4cOSpKioKAUHB5dpLyQkxOG+N2/erAkTJshqteqOO+5Q7dq1tXv3bkVHR2vbtm1auHAhgQ0AADikRoa106dP64033pC3t7cWL16sqKgoSVJBQYEmT56sr776SjNnztTChQvt95SEtddee0233HLLNfd96dIlTZo0SZK0aNEide3aVZKUmZmpMWPGaOvWrVq2bJlGjhx5zX0AAAD3USM3GGzfvl2FhYXq0qWLPahJkre3tyZOnChJ2rlzp/16YWGhkpOTVatWLYWFhTnV94oVK3T58mX16dPHHtQkKTAwULNmzZJ0ZdavuLjYqX4AAIB7qJFhzdPTU5KUnp5epuzChQuSpJtuusl+7ejRoyooKFCbNm1kMpmc6nvTpk2SpJ49e5Ypi4iIUPPmzXX+/HkdPHjQqX4AAIB7qJGvQTt16iSz2azDhw9r2rRpGjNmjIKCgrR//377JoIxY8bY65e8Aq1fv77mzJmjzZs368yZM7rpppvUvXt3jR8/XvXr13eo759++kmS1Lx5+QsSb731ViUnJysxMVHt27d35jFLsdlszNYBAPA7YbPZHK5bI8NagwYNNG/ePL300kuKjY1VbGysveymm27S/Pnz9cc//tF+rSSsffrppwoICNAdd9yhkJAQHT58WLGxsdqwYYOWLFmiFi1aXLXfS5cuKS8vT5IqDHclmxfKm/VzRm5urvbv3+/SNgEAQPWrka9BJal169bq2bOnvLy81L59e/3P//yPgoODdfHiRX3wwQdKS0uz1z1y5IgkqVevXvruu++0cOFCLVmyRHFxcerVq5cyMjL0zDPPqKio6Kp95uTkSLqyNs7Do/x/tb6+vqXqAgAAXE2NnFlLSEjQ8OHD5ePjo48//litW7eWdGUjwZtvvql//OMfGjp0qL788ktZLBZ99NFHSktLU9OmTeXt7W1vp1atWpo1a5Z+/PFHnThxQt9//73uvffeCvstCWiOrHuryvSnI/z8/Cqd+QMAAMaQlJSk3Nxch+rWyLA2Y8YMXbx4UQsWLLAHNUkym82aPHmyDhw4oL179+qTTz7R448/Ll9fXzVr1qzctvz9/XXXXXfp888/V3x8/FXDmr+/vyQpPz9fVqu13Nm1ktekfn5+zjxiGSaTyb6xAgAAGFtVNjTWuNeg+fn52rdvn0wmk7p06VKm3GQyqVu3bpKkQ4cOOdRmyYG4lSXggIAABQQESJLOnz9fbp2StWrlHbwLAADwWzUurGVlZclqtV51pqnkelFRkc6dO6cpU6bomWeeqXA35dmzZyVd2bhQmZJdoEePHi23vOQ6rywBAIAjalxYq1u3rgIDA2W1WvXtt9+WW2fr1q2SpFatWqlWrVpat26d1q9fr127dpWpm5mZqc2bN8tkMumee+6ptP+SWbv169eXKUtJSVFycrKCgoJ02223VeGpAACAu6pxYc3Dw0MDBw6UJM2cOVPJycn2MqvVqvnz52vbtm2qXbu2+vfvL4vFoj59+kiSXn31VZ0+fdpe/9KlS3r22WeVnZ2thx56SOHh4fay7OxspaSkKDU1tVT//fr1U0BAgFavXq24uDj79czMTL344ouSpBEjRvBtUAAA4BCTzdXbEg2gsLBQTz/9tDZv3iwPDw916NBBderUUWJiok6fPi2LxaIFCxaoc+fOkqTLly9r+PDhOnjwoHx9fdWhQwf5+vpq9+7dys7O1u23365FixbZNxBI0po1azRlyhSFhobav1pQ4osvvtDkyZNls9nUvn171a1bV7t371ZmZqbuvvtuvffeey4LawkJCcrJyZHFYlGrVq1c0iYAALi+qvL3u0aGNenK0Rhr1qzRmjVrlJiYqPz8fAUHB6tLly4aPXq0mjRpUqp+QUGBPvroI61du1bHjx+Xp6enbrnlFj300EMaPHhwmXB1tbAmSXv27NHChQu1f/9+FRUVqXHjxurXr58GDRpU6ngQZ+3fv1/FxcUymUwu32EKAACuj9zcXNlsNnl6eqpdu3ZXrVtjw5q72Ldvn8vPbAMAADeGyWRShw4drlqnRp6z5k7MZrMKCwvl4eEhHx+f6h4OAABwQMmZrI4si2JmDQAAwMBq3G5QAACAmoSwBgAAYGCENQAAAAMjrAEAABgYYQ0AAMDACGsAAAAGRlgDAAAwMMIaAACAgRHWAAAADIywBgAAYGCENQAAAAMjrAEAABgYYQ0AAMDACGsAAAAGRlgDAAAwMMKam3vxxRfVokULRUVF6cKFC1etW1RUpH79+qlFixYaP378DRohAADX386dO9WiRQu1aNFCGzduvGrdNWvWqEWLFho7duwNGRthzc29+OKLCg0NVWZmpqZNm3bVuu+//74OHz6sm2++WTNmzLhBIwQA4MZ6+eWXlZGRUd3DsCOsubmAgADNnj1bHh4e2rhxo7744oty6yUkJOj999+XyWTS7NmzFRQUdINHCgDAjfHzzz/r5Zdfru5h2BHWoDvvvFPDhg2TJM2YMaPM69CCggL97W9/U2FhoR5//HF17dq1OoYJAMB1FxgYKF9fX8XFxWn16tXVPRxJhDX8f88995xuvfXWcl+HLliwQElJSWrevLkmTZpUTSMEAOD6u/nmm/X8889Lkv7+97/r9OnT1Twiwhr+Px8fH82ZM0dms1kbN260L65MTEzU4sWL5ePjozfeeEPe3t7VPFIAAK6voUOHqlOnTrp8+bJeeOEF2Wy2ah0PYQ12bdq00ZNPPinpyuvQy5cva+rUqSoqKtKkSZPUvHnzah4hAADXn8lk0qxZs1SrVi3t2rVLS5curdbxENZQypNPPqnIyEidPXtWAwcOVHx8vO6++24NHTq0uocGAMANExISoqlTp0qS5s2bp6NHj1bbWAhrKMXLy0tz5syRr6+vkpOTFRQUpNmzZ8tkMlX30AAAuKH69Omjnj17Kj8/X5MnT1ZhYWG1jIOwhjIiIiL0xz/+UZL06KOPql69etU8IgAAqsdrr72mevXq6fDhw1qwYEG1jIGwhnJ5eXmV+l8AANxRUFCQpk+fLkn64IMPdODAgRs+BsIaAADAVdx7773q37+/iouLNXnyZOXm5t7Q/glrAAAAlSj5POOJEyf0wQcf3NC+CWsAAACVCAgI0Jw5c+Th4aH//Oc/N7RvwhoAAIADOnbsaP88441EWAMAAHDQ//7v/6pZs2Y3tE+Trbq/oQAAAIAKMbMGAABgYIQ1AAAAAyOsAQAAGBhhDQAAwMAIawAAAAZGWAMAADAwwhoAAICBEdYAAAAMjLAGAABgYIQ1AKhmubm5OnXq1DXd+9NPP7l4NACMhrAGANXoiy++UK9evbR9+/Yq3Xf8+HGNHDlSr7zyynUaGQCjIKwBQDWaN2+ezp07V+X71q5dqy1btlyHEQEwGsIaAACAgRHWAAAADMxks9ls1T0IAHA37777rubPn1/m+oQJExQaGqopU6boz3/+s4YMGaLXX39dKSkpCgwM1KhRozRr1qwy94WGhmrTpk03YugAbjBm1gCgGoSEhKhDhw7y9vaWJDVt2lQdOnRQSEiIvc6xY8c0atQonT59Ws2aNVNWVpYCAwNL1QsICFCHDh3Upk2bankOANcfM2sAUI3uvfdenT59WjNmzNCAAQMkSWvWrNGUKVMkSe3atdPixYtVq1YtXbx4UYGBgTKZTPaZuQ4dOmjVqlXV+QgArjOv6h4AAKBiEydOVK1atSRJN910UzWPBkB14DUoABiUh4eH2rdvX93DAFDNCGsAYFC1a9eWr69vdQ8DQDUjrAGAQfn4+FT3EAAYAGENAADAwAhrAAAABkZYA4BqZDKZJElVPUXpWu8D8PtDWAOAamSxWCRJp0+frtJ9/v7+kqT09HQVFRW5fFwAjIOwBgDVqHXr1pKkxYsXq2/fvoqJiXHovlatWkm6EvJ69uypxx57jFk2oIYirAFANfrb3/6mXr16yc/PT8eOHVNKSopD9911112aPHmyQkNDlZ6erlOnTunChQvXebQAqgOfmwIAADAwZtYAAAAMjLAGAABgYIQ1AAAAAyOsAQAAGBhhDQAAwMAIawAAAAZGWAMAADAwwhoAAICBEdYAAAAMjLAGAABgYIQ1AAAAAyOsAQAAGBhhDQAAwMD+H48lKCdLnR5YAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Plot\n", "ax = sns.pointplot(data=p,\n", " x='trt', y='estimate',\n", " linestyle='none')\n", "# Add errors\n", "ax.vlines(p['trt'], p['conf_low'], p['conf_high'])" ] }, { "cell_type": "markdown", "id": "fb0f060e-cf5a-48d2-adb1-11adace10533", "metadata": {}, "source": [ "## Testing hypotheses" ] }, { "cell_type": "markdown", "id": "6d2c0b0b-ffaa-44ce-a4ad-c7ad12bd46af", "metadata": {}, "source": [ "We can compute the difference between these two as usual with the `hypothesis` argument, and obtain the interval around the difference:" ] }, { "cell_type": "code", "execution_count": 6, "id": "65c3f5ef-9e4a-4569-a1b1-2a76eff9bb5b", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "shape: (1, 8)
termestimatestd_errorstatisticp_values_valueconf_lowconf_high
strf64f64f64f64f64f64f64
"b1=b2"9.4902941.496986.3396282.3032e-1032.0156416.55626812.42432
" ], "text/plain": [ "shape: (1, 8)\n", "┌───────┬──────────┬───────────┬──────┬─────────┬─────┬──────┬───────┐\n", "│ Term ┆ Estimate ┆ Std.Error ┆ z ┆ P(>|z|) ┆ S ┆ 2.5% ┆ 97.5% │\n", "│ --- ┆ --- ┆ --- ┆ --- ┆ --- ┆ --- ┆ --- ┆ --- │\n", "│ str ┆ str ┆ str ┆ str ┆ str ┆ str ┆ str ┆ str │\n", "╞═══════╪══════════╪═══════════╪══════╪═════════╪═════╪══════╪═══════╡\n", "│ b1=b2 ┆ 9.49 ┆ 1.5 ┆ 6.34 ┆ 2.3e-10 ┆ 32 ┆ 6.56 ┆ 12.4 │\n", "└───────┴──────────┴───────────┴──────┴─────────┴─────┴──────┴───────┘\n", "\n", "Columns: term, estimate, std_error, statistic, p_value, s_value, conf_low, conf_high" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Difference\n", "p = me.predictions(cog_mod,\n", " newdata=me.datagrid(cog_mod,\n", " trt=['Y', 'N'],\n", " age=2),\n", " hypothesis='b1=b2')\n", "p" ] }, { "cell_type": "code", "execution_count": 7, "id": "a6fba33d-7132-430a-b912-b2b3f3f9011d", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAksAAAHHCAYAAACvJxw8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAuj0lEQVR4nO3de5iN9f7/8de9ZoaZwdCMZjBOsTPIyBRKir5KUruQaEeHbdux00GuDspX+Trs3W475kwqRqKTL5U2IpUioS8Nxsw2gxjtcRyGWczMWvfvj37WbjI+Zqw1a1nL83FdXVdr3Z+573eua+U597rXvSzbtm0BAACgVI5ADwAAAHApI5YAAAAMiCUAAAADYgkAAMCAWAIAADAglgAAAAyIJQAAAIPwQA8Q7NLS0lRUVCSHw6HKlSsHehwAAFAGZ86ckdvtVkREhJKTk41riSUvFRUVybZtuVwuFRQUBHocAABQDkVFRRdcQyx5yeFwyOVyybIsRUVFBXocAABQBk6nU7Zty+G48BVJxJKXKleurIKCAkVFRalZs2aBHgcAAJRBenq6CgoKynQJDRd4AwAAGBBLAAAABsQSAACAAbEEAABgQCwBAAAYEEsAAAAGxBIAAIABsQQAAGBALAEAABgQSwAAAAbEEgAAgAGxBAAAYEAsAQAAGBBLAAAABuGBHgAALkWni1ya8WWW5/HjtzZWZERYACcCECjEEgCU4kyxW6+v/pfn8Z9uvopYAi5TvA0HAABgQCwBAAAYEEsAAAAGxBIAAIABsQQAAGBALAEAABgQSwAAAAbEEgAAgAGxBAAAYEAsAQAAGBBLAAAABsQSAACAAbEEAABgQCwBAAAYEEsAAAAGQRlLe/bsUatWrTRq1KhStzudTk2fPl3dunVTq1at1LJlS3Xt2lXjxo1TXl6ef4cFAABBLTzQA5TX4cOHNWjQIDmdzlK35+Xl6eGHH1ZmZqaqVaumlJQUhYWFKS0tTW+88YY+++wzvfvuu6pVq5afJwcAAMEoqM4spaenq0+fPsrKyjrvmrFjxyozM1Nt27bVypUr9fbbb2vOnDn6/PPPdcsttygnJ0cvv/yyH6cGAADBLChi6fjx4xo7dqx69+6tvXv3qm7duqWuO336tD799FNJ0t///nfFxsZ6tsXExOi1116TZVlau3atjh075pfZAQBAcAuKWEpNTdWcOXMUGxurGTNmqHv37qWuO3LkiK655hpdd911SkxMPGd7XFycqlevLtu2dfDgwQqeGgAAhIKguGapVq1aGjp0qPr06aPIyEht37691HWJiYl69913z7ufvXv3Ki8vTw6HQwkJCRU1LgAACCFBEUu9evXyyX7GjRsnSbrppptUo0YNn+zzLNu25XK5fLpPAIHj/s3r2e1yyeUKipPxAMrAtu0yrw2KWPKFmTNnauXKlYqMjNQLL7zg8/07nU5t2bLF5/sFEBinCt0lHqelpalKJWIJuBxdFrE0efJkTZs2TQ6HQ3/729+UlJQU6JEAAECQCOlYKiws1PDhw7V06VKFh4fr1Vdf1d13310hx4qKiiLCgBBywlkkLV3teZycnKyYqIgATgTAlzIyMs57z8bfCtlYOnLkiJ544gn93//9n6pWrarXX39dN998c4Udz7IshYWFVdj+AfiXI8z9m8dhvMaBEGJZVpnXhmQs/fTTT/rjH/+onJwc1a1bVzNnztTVV18d6LEAAEAQCrlYys3N1SOPPKKff/5ZLVu21MyZMxUXFxfosQAAQJAKuVh6/vnn9fPPPyspKUlz585VlSpVAj0SAAAIYiEVS99++602bNggSapWrZpGjBhx3rWDBw9WvXr1/DUaAAAIUiEVS2vWrPH8+6ZNm4xrH330UWIJAABckGWX5xaWOEd6eroKCgoUHR2tZs2aBXocAD5y3Fmka0eu9DzeOuIOVefWAUDIKM/f39yOFgAAwIBYAgAAMCCWAAAADIglAAAAA2IJAADAgFgCAAAwIJYAAAAMiCUAAAADYgkAAMCAWAIAADAglgAAAAyIJQAAAANiCQAAwIBYAgAAMCCWAAAADIglAAAAA2IJAADAgFgCAAAwIJYAAAAMiCUAAAADYgkAAMCAWAIAADAglgAAAAyIJQAAAANiCQAAwIBYAgAAMCCWAAAADIglAAAAA2IJAADAgFgCAAAwIJYAAAAMiCUAAAADYgkAAMCAWAIAADAglgAAAAyIJQAAAANiCQAAwIBYAgAAMCCWAAAADIglAAAAg6COpT179qhVq1YaNWqUT9YBAAD8VtDG0uHDhzVo0CA5nU6frAMAAChNUMZSenq6+vTpo6ysLJ+sAwAAOJ/wQA9QHsePH9fs2bOVmpqqwsJC1a1bV/v377/odQAAABcSVGeWUlNTNWfOHMXGxmrGjBnq3r27V+sAAAAuJKhiqVatWho6dKhWrFihTp06eb0OAADgQoLqbbhevXr5dJ0v2bYtl8vl9+MCqBju37ye3S6XXK6g+v0SgIFt22VeG1SxdClzOp3asmVLoMcA4COnCt0lHqelpalKJWIJuBzxygcAADDgzJKPREVFKSkpKdBjAPCRE84iaelqz+Pk5GTFREUEcCIAvpSRkVHmezASSz5iWZbCwsICPQYAH3GEuX/zOIzXOBBCLMsq81rehgMAADAglgAAAAyIJQAAAANiCQAAwIBYAgAAMLDs8tzCEudIT09XQUGBoqOj1axZs0CPA8BHjjuLdO3IlZ7HW0fcoercOgAIGeX5+5szSwAAAAbEEgAAgAGxBAAAYEAsAQAAGBBLAAAABsQSAACAAbEEAABgQCwBAAAYEEsAAAAGxBIAAIABsQQAAGBALAEAABgQSwAAAAbEEgAAgAGxBAAAYEAsAQAAGBBLAAAABsQSAACAAbEEAABgQCwBAAAYEEsAAAAGxBIAAIABsQQAAGBALAEAABgQSwAAAAbEEgAAgAGxBAAAYEAsAQAAGBBLAAAABsQSAACAAbEEAABgQCwBAAAYEEsAAAAG4b7a0fHjx7V+/XplZ2crPz9fQ4cO1ZkzZ7R161a1bdvWV4cBAADwK69jybZtTZkyRW+//bZOnz7teX7o0KHav3+/Hn30UbVq1UrTpk1TbGyst4cDAADwK6/fhnvhhRc0Y8YMOZ1O1ahRQ1FRUZ5teXl5sm1bW7Zs0cMPPyyn0+nt4QAAAPzKq1hauXKlPvnkE8XGxuqNN97Q+vXr1bRpU8/266+/XgsWLFBcXJyys7OVmprq9cAAAAD+5FUsvffee7IsS+PHj9ctt9xS6prrr79ekyZNkm3bWrFihTeHAwAA8DuvYmnbtm2qXbu2brzxRuO61q1bKzExUXv27PHmcAAAAH7nVSwVFBSoRo0aZVobGxur4uJibw4HAADgd17FUs2aNbV3717Ztm1cV1RUpD179qhmzZreHA4AAMDvvIqltm3bqqCgQO+++65x3bx585Sfn6/WrVt7cziPPXv2qFWrVho1atR516xbt079+vVTu3btlJKSop49e+qDDz64YNgBAAD8mlex1K9fPzkcDr322mtKTU3VsWPHSmw/cuSIJk2apAkTJsjhcOihhx7yalhJOnz4sAYNGmS8DcHChQvVr18/bdy4Uc2bN9cNN9ygrKwsDR8+XMOGDfN6BgAAcPnw6qaUTZs21bBhwzRmzBi9+uqrevXVVz3b2rVrp7y8PEm/3Lhy8ODBatmypVfDpqena/Dgwdq7d+951+zevVujR49W1apVNX/+fDVv3lySdODAAT366KNavHixOnbsqDvvvNOrWQAAwOXB65tS9u3bV7NmzVJSUpJs2/b8c+zYMdm2rfr162vChAl6/PHHL/oYx48f19ixY9W7d2/t3btXdevWPe/aOXPmyOVyqX///p5QkqQ6derolVde8awBAAAoC598N1yHDh3UoUMH5eTk6F//+pfy8/MVFRWlq666So0bN/Z6/6mpqZozZ45q1aqlESNGaPv27Zo6dWqpa9esWSNJuuOOO87ZdtNNN6latWpKS0tTbm6uEhISvJ4NAACENq9iaePGjapWrZrnrt2JiYlKTEwsde0333yjvXv3qm/fvuU+Tq1atTR06FD16dNHkZGR2r59e6nrDh8+rCNHjigiIkKNGjU6Z3tYWJgaNWqkrVu3KiMjw6exZNu2XC6Xz/YHILDcv3k9u10uuVxen4wHcIkozwe+vIqlhx9+WK1bt9Y777xzwbWTJk266Fjq1atXmdYdPHhQ0i+3NHA4Sv+fWnx8fIm1vuJ0OrVlyxaf7hNA4JwqdJd4nJaWpiqViCXgclTmWDp58uQ5n3aTpNOnT2vfvn3n/TnbtpWTk6Ps7OwK/9h+QUGBJCkyMvK8aypXrlxiLQAAgEmZY+nUqVO65557dObMGc9zlmVp+/btpV4fVJpWrVqVe8DyOHs2ybKsC671dbhFRUUpKSnJp/sEEDgnnEXS0tWex8nJyYqJigjgRAB8KSMjw3gbol8rcywlJCSoX79+mjFjhuc5y7LKHB116tTR8OHDy3q4i1KlShVJv5ztOp+zsRcdHe3TY1uWpbCwMJ/uE0DgOMLcv3kcxmscCCFlObFyVrmuWRo0aJDuv/9+Sb+cmbn99tuVnJysSZMmnfdnHA6HoqOjVb169fIc6qKcvWD78OHDsm271D+Is9cqnb12CQAAwKRcsRQREVHi025t2rRRUlLSeT8B5281atRQQkKCcnNztXfvXjVs2LDEdpfLpezsbEniLTMAAFAmXn20Y/78+RX+1lp5dezYUZK0cuXKc7Z9++23ys/PV9OmTVWrVi1/jwYAAIKQzz8H63a7S/xTWFioEydOKDMzU7NmzfL14c7Rt29fhYWFadasWSU+yn/gwAGNHj1akjRw4MAKnwMAAIQGr+/g/fXXX2vKlCnKzMxUYWHhBddXdKg0bdpUzzzzjMaPH68+ffqobdu2ioyM1IYNG1RQUKBevXrprrvuqtAZAABA6PAqlrZt26ZBgwbJ5XJd8FNx4eHhSklJ8eZwZTZgwAA1btxYc+fO1Y8//ijLstS4cWM9+OCD6tGjh19mAAAAocGrWEpNTVVxcbGuvvpq/fnPf1ZkZKQGDx6sLl266IEHHtC///1vLV68WJs2bVLr1q01d+5cnwz91FNP6amnnjKuue2223Tbbbf55HgAAODy5VUsbdq0SWFhYZo8ebKuuuoqSVLt2rW1b98+3XTTTZKkHj166Omnn9aqVau0bNky3X333d5PDQAA4CdeXeB95MgR1alTxxNK0i/XDP36+iXLsvTSSy9JkhYvXuzN4QAAAPzO60/D1ahRo8Tjhg0byuVyaffu3Z7n6tSpowYNGigjI8PbwwEAAPiVV7EUFxenQ4cOlXiubt26kqRdu3aVeL5KlSrKy8vz5nAAAAB+51UstWjRQrm5uVq7dq3nuUaNGsm2bW3cuNHz3OnTp/XTTz/55StPAAAAfMmrWOrWrZts29bTTz+tsWPHqri4WCkpKapevbo+/PBDLVmyRJmZmRo+fLjy8/NLXNsEAAAQDLyKpdtuu0133323nE6n5s2bp7CwMEVGRurBBx9UcXGxXnrpJXXr1k3Lli2TZVl65JFHfDU3AACAX3h9B+/x48erffv2+vbbb2VZliTpySefVG5urpYsWSLbthUWFqZ+/frpjjvu8HpgAAAAf/I6liTpvvvu03333fefnYaH69VXX9WQIUN04MAB1a9fX7Gxsb44FAAAgF/5JJbOJz4+XvHx8RV5CAAAgArlk1gqLCxUdna2Tp48ecHviGvTpo0vDgkAFebgidN6f9O+Es/NX79HvVvXU3xMZICmAhAoXsfSxIkTlZqaqtOnT19wrWVZ2rFjh7eHBIAKkbb/uGZ+naUV2/6tYnfJX/zGrczUpFX/UpcWtfSXDo2VXJdboQCXC69iae7cuZo1a1aZ11/orBMABMo73+3VK0u3yW3431Sx29ayH3/WP9N+1qhuLfTQjQ38NyCAgPEqlt5//31ZlqX77rtPTzzxhOLj4xUeXqGXQQGAz73z3V4NX7KtzOvdtjR8yTZZltT3BoIJCHVe3Wdp3759iouL0+jRo1WnTh1CCUDQSdt/XK8sLXso/drLS7Ypbf9xH08E4FLjVSxVrVpVV155pRwOr7+PFwACYubXWca33kzctjTr6yzfDgTgkuNV5bRp00a7d+/WqVOnfDUPAPjNwROntWLbv73ax/Jt/9bBExf+gAuA4OVVLA0aNEgul0sjR47k4m0AQefjrQfO+dRbeRW7bX289YCPJgJwKfLqIqOmTZvq9ddf15NPPqktW7aoffv2io2N9XztSWmefPJJbw4JAD6Tdeikj/bD2XUglHkVSydPntTs2bPldru1b98+LVq06II/QywBuFScOuPy0X6KfbIfAJcmr2Jp4sSJ2rJliyTpyiuvVJ06dRQREeGLuQCgwlWpHOaj/fBJYCCUefUKX716tSzL0t/+9jf16NHDVzMBgF80vrKqj/ZTxSf7AXBp8uoC76NHj6pBgwaEEoCgdO+1dRTuOP81lmUR7rB077V1fDQRgEuRV7EUHx+vsDDfnMYGAH+Lj4lUlxa1vNrHnS1q8eW6QIjzKpa6dOmi7Oxsbd++3VfzAIBf/aVDY13sySWHJQ3s0Ni3AwG45HgVS48//rjq16+vxx9/XMuXL+fmlACCTnLd6hrVrcVF/ezo7i2UXLe6jycCcKnx6gLvUaNGqV69evrmm280ZMgQWZalatWqKSoqqtT1lmVpzZo13hwSAHzuoRsbyLJ++a63styj0mH9Ekp8iS5wefAqlj7++GPPv9u2Ldu2dfz4cR0/XvoXS5puVgkAgdT3hgZqmVhDs77O0vJt/y71zt7hDkt3tqilgR0ac0YJuIx4FUuvvvqqr+YAgIBLrltdU/tcp4MnTuv9Tfs0bmWmZ9tzdzRR7zb1FF+Ni7mBy41XscQtAwCEoviYSD3crmGJWHq4XUNVj+Kmu8DlyKsLvAEAAEJdmc8srV+/XpJ03XXXqXLlyiWeK4927dqV+2cAAAACpcyx1K9fPzkcDi1btkxXXXWV57nyXLRtWZZ27NhR/ikBAAACpFzXLLnd7nOes+0yfM72ItYCAABcCsocSzt37izTcwAAAKHErxd4HzlyxJ+HAwAA8JpXsXTbbbdpyJAhZVr7hz/8gVsNAACAoONVLOXk5OjgwYMXXOd2u3Xo0CEdO3bMm8MBAAD4XZmvWdq1a5dGjBhxzvOZmZnq27fveX/Otm3l5ubqwIEDqlOnzsVNCQAAECBljqXf/e53ioyM1Lfffut5zrIs5efna/PmzWXax0MPPVT+CQEAAAKoXLcOePnll/Xpp596Hk+dOlV16tTRfffdd96fsSxLVapUUbNmzXTDDTdc/KTl9M9//lPvvvuuduzYIdu21bBhQ/Xs2VO9evVSpUqV/DYHAAAIbuWKpYYNG+rJJ5/0PJ46dapq165d4rlLwYgRI7Ro0SJJUpMmTVS3bl3t3LlTo0aN0meffabp06erenW+MRwAAFyYV1+ku3r1as9Xn1wqli5dqkWLFikiIkJjx45V165dJUlFRUV67bXXNH/+fI0ZM0Zjx44N8KQAACAYePVpuMTERNWsWbPUbadPn9YXX3yhVatWKS8vz5vDlMvChQslSY899pgnlCQpIiJCL774oho3bqyPP/5YmZmZ59sFAACAh1dnliQpNzdXM2bMUJ06dTRgwABJUlZWlvr166dDhw5JkqKiojRmzBjddddd3h7ugjIyMiRJt99++znbwsPD1aZNG2VlZemrr75SkyZNKnweAAAQ3LyKpaNHj6p37946ePCgbr31Vs/zr7zyig4ePOi5uPvkyZN64YUXlJSUpMaNG3s7s5HL5ZIkVa1atdTt4eG//CdnZWX59Li2bXuODSD4uX/zena7XHK5/PqlBwAqUHm+r9arWJo3b55yc3PVoEEDPfDAA5KkvXv3avPmzQoLC9OCBQvUqlUrTZgwQbNnz9bcuXM1evRobw55QY0aNVJ6erq+//57NWjQoMQ227b1ww8/SPL9V684nU5t2bLFp/sEEDinCkt+cXhaWpqqVCKWgMuRV6/8r7/+WuHh4XrzzTc9Z5a+/PJLSdJ1112nVq1aSZKeeuopxcTE6LvvvvPmcGVy9jYG48aNKxEvbrdbkydP1o4dOyRJhYWFFT4LAAAIfl6dWdq3b58aNmyounXrep5bt26dLMvSTTfd5HkuIiJCdevW9flbX6Xp27evvvvuO61evVoPPvigkpOTVbNmTWVkZCg3N1d/+MMftGjRIs/bcb4SFRWlpKQkn+4TQOCccBZJS1d7HicnJysmKiKAEwHwpYyMDDmdzjKt9aoYXC5XiRs8FhcXa+PGjZKktm3blljrdDplWZY3hyuTsLAwTZ06Ve+++64++OADpaenKzo6WjfccIOmTp2q3bt3a9GiRT6/z5JlWQoLC/PpPgEEjiPM/ZvHYbzGgRBSnibxKpYSExOVk5OjoqIiRUREaOPGjSooKFDVqlU9b8FJv3xibt++fapXr543hyszh8Ohhx56qNSvV/n88889swMAAFyIV9csJScn68SJExo3bpx27typSZMmybIsdezY0fMb2JEjR/T888/L5XKpXbt2Phna5KefftI333zjuW3Bb61bt06S1LJlywqfBQAABD+vYumxxx5TZGSkUlNT1aNHD23dulVhYWF67LHHJEmbNm1Sx44dtXHjRlWrVk1/+tOffDK0yUcffaT+/fvrww8/PGfbjh07tGXLFtWoUUPt27ev8FkAAEDw8yqWGjVqpLfeekvJycmqVKmSmjRpohkzZqhp06aSpPj4eBUXF+vqq6/WwoULS1wIXlFuv/12WZaluXPnat++fZ7nf/75Zz333HOybVsDBw5UdHR0hc8CAACCn9cfCUtJSdH7779f6ra6detqyZIlnnjyh+TkZP35z3/WG2+8oXvuuUdt2rSRJG3YsEFnzpxRjx499Mc//tFv8wAAgODms8/Pr127VmvWrFF2drby8/P10Ucf6eTJk1q1apUSEhJ0xRVX+OpQF/Tss8+qXr16Wrhwob777jtVqVJFLVu2VJ8+fdS1a1e/fCoPAACEBq9j6ciRI3rmmWe0adMmSb/cJftsjBw4cEBTp07V/PnzNXv2bF177bXeHq5MLMvSAw884LmrOAAAwMXy6pqlwsJC9e/fXxs3blSVKlXUuXNnJSQk/GfnDodq1Kih48ePq1+/fsrJyfF6YAAAAH/yKpYWLFignTt3qlWrVlq5cqUmT55c4v5FTZo00apVq5SSkiKn06m3337b64EBAAD8yatYWrZsmRwOh8aOHavY2NhS11StWlXjxo1TWFiY1q5d683hAAAA/M6rWMrOzlbjxo0veGfuxMRENWzYUD///LM3hwMAAPA7r2LJ7XZfeNH/FxERwfcqAQCAoONVLCUmJmrPnj06efKkcd2xY8f0r3/9i+9jAwAAQcerWOrYsaOKioo0duxY47oxY8bI5XLplltu8eZwAAAAfufVfZb69++vjz76SO+//76OHDmie+65R/n5+ZKkrKwsZWZmasGCBdq8ebOqVKnCnbMBAEDQ8SqW4uLiNH36dA0aNEirVq3S6tWrPdt+//vfS/rlJpXR0dGaMGFCiXswAQAABAOv3oaTpOuvv14ff/yxHnnkEdWuXVu2bXv+iYuL0/33368lS5aoQ4cOvpgXAADAr3zy3XAJCQkaNmyYhg0bpoKCAuXn5ys6OlrVqlXzxe4BAAACxmdfpHtWdHS0oqOjfb1bAACAgPD6bTgAAIBQRiwBAAAYEEsAAAAGxBIAAIABsQQAAGBALAEAABgQSwAAAAbEEgAAgAGxBAAAYEAsAQAAGBBLAAAABsQSAACAAbEEAABgQCwBAAAYEEsAAAAGxBIAAIABsQQAAGBALAEAABgQSwAAAAbEEgAAgAGxBAAAYEAsAQAAGBBLAAAABsQSAACAAbEEAABgQCwBAAAYEEsAAAAGxBIAAIBBeKAHqEhfffWV5s6dq7S0NJ0+fVoJCQm69dZb9fjjj6tmzZqBHg8AAASBkD2z9Pbbb2vAgAFav369rr76anXs2FGFhYV655131L17d/3000+BHhEAAASBkIylnJwcjR8/XpUqVdK8efO0cOFCTZs2TatXr1bXrl116NAh/fWvfw30mAAAIAiEZCytX79eRUVFat++vW644QbP85UqVdIzzzwjSdqwYUOApgMAAMEkJGMpLCxMknTw4MFzth0+fFiSdMUVV/h1JgAAEJxCMpbatWuniIgIbd++XSNGjFBOTo6cTqfWr1+vl156SZI0YMCAAE8JAACCQUh+Gq5WrVqaOHGi/vu//1uLFi3SokWLPNuuuOIKTZ06VZ07d/bpMW3blsvl8uk+AQSO+zevZ7fLJZcrJH+/BC5Ltm2XeW1IxpIkNW/eXHfccYf+93//V8nJyapRo4a2b9+ugwcPavbs2WratKnq1avns+M5nU5t2bLFZ/sDEFinCt0lHqelpalKJWIJuByFZCylp6erX79+qly5sj744AM1b95cklRUVKQJEyborbfe0sMPP6zPPvtM0dHRAZ4WAABcykIylsaMGaNjx45p2rRpnlCSpIiICL3wwgvaunWrNm/erA8//FCPPPKIT44ZFRWlpKQkn+wLQOCdcBZJS1d7HicnJysmKiKAEwHwpYyMDDmdzjKtDblYOnPmjH744QdZlqX27dufs92yLHXs2FGbN2/Wtm3bfHZcy7I8n8IDEPwcYe7fPA7jNQ6EEMuyyrw25N6AP3HihNxutzFezj5fXFzsz9EAAEAQCrlYiouLU40aNeR2u/Xll1+Wuubbb7+VJDVr1syPkwEAgGAUcrHkcDj04IMPSpL++te/KjMz07PN7XZr6tSpWrdunWJiYtSzZ89AjQkAAIJEyF2zJElPPPGEdu7cqTVr1qhbt2667rrrVL16de3cuVM5OTmKjo7W66+/rtjY2ECPCgAALnEhGUsRERGaMWOGFi9erMWLF2vnzp06c+aM4uPj1bt3bz322GOqX79+oMcEAABBICRjSfrlKveePXvyVhsAAPBKyF2zBAAA4EvEEgAAgAGxBAAAYEAsAQAAGBBLAAAABsQSAACAAbEEAABgQCwBAAAYEEsAAAAGxBIAAIABsQQAAGBALAEAABgQSwAAAAbEEgAAgAGxBAAAYEAsAQAAGBBLAAAABsQSAACAAbEEAABgQCwBAAAYEEsAAAAGxBIAAIABsQQAAGBALAEAABgQSwAAAAbEEgAAgAGxBAAAYEAsAQAAGBBLAAAABsQSAACAAbEEAABgQCwBAAAYEEsAAAAGxBIAAIABsQQAAGBALAEAABgQSwAAAAbEEgAAgAGxBAAAYEAsAQAAGIQHegBf69Spk3Jyci64rm3btpo/f74fJgIAAMEs5GLp9ttv19GjR0vdZtu2li9fruLiYl1zzTV+ngwAAASjkIulYcOGnXfbtGnTVFxcrDZt2ui5557z41QAACBYXTbXLG3YsEFTp05VTEyMxo8fr/DwkOtEAABQAS6LWCosLNQrr7wit9utF154QQkJCYEeCQAABInLIpbmzp2rPXv2KDk5Wffff3+gxwEAAEEk5N+LOnnypGbPni1Jevrpp2VZVoUcx7ZtuVyuCtk3AP9z/+b17Ha55HJdFr9fApcF27bLvDbkY2nRokXKz8/XNddcow4dOlTYcZxOp7Zs2VJh+wfgX6cK3SUep6WlqUolYgm4HIX0K9/lcik1NVWSNGDAgABPAwAAglFIn1n6/vvvlZubq+rVq6tTp04VeqyoqCglJSVV6DEA+M8JZ5G0dLXncXJysmKiIgI4EQBfysjIkNPpLNPakI6lFStWSJK6dOmiSpUqVeixLMtSWFhYhR4DgP84wty/eRzGaxwIIeW5hjmk34b76quvJEl33nlngCcBAADBKmRj6dChQzpw4IDCw8OVkpIS6HEAAECQCtlY+vHHHyVJTZo0UXR0dICnAQAAwSpkY2nfvn2SpHr16gV4EgAAEMxCNpaOHTsmSapdu3aAJwEAAMEsZD8NN2TIEA0ZMiTQYwAAgCAXsmeWAAAAfIFYAgAAMCCWAAAADIglAAAAA2IJAADAgFgCAAAwIJYAAAAMiCUAAAADYgkAAMCAWAIAADAglgAAAAyIJQAAAANiCQAAwIBYAgAAMCCWAAAADIglAAAAg/BADwAAl6LK4Q4Nvu3qEo8BXJ6IJQAoRWREmIZ0bhLoMQBcAvhVCQAAwIBYAgAAMCCWAAAADIglAAAAA2IJAADAgFgCAAAwIJYAAAAMiCUAAAADYgkAAMCAWAIAADAglgAAAAyIJQAAAANiCQAAwIBYAgAAMAgP9ADB7syZM5Ikp9Op9PT0AE8DAADKwul0SvrP3+MmxJKX3G63JMm2bRUUFAR4GgAAUB5n/x43IZa8FBERoaKiIjkcDlWuXDnQ4wAAgDI4c+aM3G63IiIiLrjWsm3b9sNMAAAAQYkLvAEAAAyIJQAAAANiCQAAwIBYAgAAMCCWAAAADIglAAAAA2IJAADAgFgCAAAwIJYAAAAMiCUAAAADYgkAAMCAWAIAADAglgAAAAyIJQAAAANiCQAAwIBYAhAUpkyZoqSkJI0aNeqi97F06VIlJSVp+fLlPpvrxRdfVFJSkt58881y/dzRo0f1j3/8Q3fddZdatmypVq1aqVu3bpo5c6ZOnz7ts/kAeC880AMAgD/88MMPGjlyZKDHkCTt27dPffv2VW5urmJjY3XjjTeqsLBQW7du1cSJE7Vy5UqlpqaqatWqgR4VgDizBOAy8Omnn6p///46depUoEeRJL388svKzc1V165dtXr1as2ePVtz587V8uXLdc0112j79u0aN25coMcE8P8RSwBC1u7du/X000/r2WeflSTVrFkzwBNJ+/fv1/r161WtWjWNGTNG0dHRnm0JCQn6n//5H0m/BJ5t2wGaEsCv8TYcgKCzbt06TZkyRTt27FBkZKTatGmjgQMHKjk5ucS6V155Rd9//71SUlI0ZswYjRw5UocPHy51nxs2bNAjjzxSpuP36NFDf//73895ftmyZZozZ4527dqlmJgYtW/fXoMGDVLDhg09a44cOaKUlBQlJCSU+jZbo0aNJEn5+fk6deoUb8UBlwBiCUBQWbdunRYuXKhatWqpY8eO2rdvnz7//HOtWbNGkyZNUufOnT1rW7RooYceekh33HGHLMsy7rdmzZq65557yjRDSkrKOc8tXrxYu3bt0lVXXaX/+q//UmZmppYuXarPP/9cb731ludnrr32Wi1atOi8+/7xxx8lSTExMYQScIkglgAEld27d6tHjx4aNWqUKlWqJElasGCBRo0apWHDhql169a64oorJElDhw4t834bN27s1XVCu3bt0l/+8hc988wzsixLtm1r/PjxeuONN/Tcc89p+fLlioiIMO7D5XJp4sSJkqSuXbte9CwAfItrlgAElbi4OL388sueUJKkvn37qmPHjjpx4oSWLl0akLmaNGmiwYMHe85gWZalZ599Vk2aNNH+/fv15ZdfGn/etm2NHDlSP/74o2JjY/Xkk0/6YWoAZUEsAQgqnTt3VpUqVc55/vbbb5ckfffdd/4eSZJ0zz33yOEo+b9Uy7LUqVMnSb9cE3U+xcXFGjZsmN577z1VrlxZkydPVnx8fIXOC6DseBsOQFCpV69eqc/Xrl1bknTw4MGL2q+3F3hfaK7c3NxSt584cUJDhgzRN998o+joaE2fPl1t2rQpx+QAKhqxBCCoVK5cudTnz37M/kLXBZ2Ptxd4X8xce/fu1cCBA7V7927Fx8dr5syZuuaaa8oxNQB/IJYABJXznaHZv3+/pP+cySkvby/wLu9cP/74ox577DHl5eWpefPmmjFjhmrVqnXRxwdQcbhmCUBQWbt2bak3a/znP/8pSbrxxhv9PZIk6euvvz7nueLiYq1atUpSybkyMjLUv39/5eXl6dZbb9WCBQsIJeASRiwBCCo7d+7UpEmTSgTTrFmz9P333yshIUH33ntvQOb64osvStw/qbi4WGPGjNGePXvUvHlz3XzzzZKkwsJCPfPMMzpx4oRuvvlmTZs2rcRdvAFcengbDkBQSUlJ0cyZM7VixQolJSVp165d2rVrl6pVq6bJkycHLDxSUlI0YsQIvffee6pfv762bdum/fv3KyEhQRMnTvTcUmDJkiXKzs6WJDkcDr344ovn3efIkSNL/eQfAP8ilgAEld///vfq16+fZs2apS+++EJVq1ZVt27d9NRTT533E2n+MHDgQOXm5mrevHlavXq1YmNj1adPHz3xxBMlvpNuzZo1nn8v7a27Xxs2bBixBFwCLJtvagQAADgvrlkCAAAwIJYAAAAMiCUAAAADYgkAAMCAWAIAADAglgAAAAyIJQAAAANiCQAAwIBYAgAAMCCWAIQMp9Op/fv3B3oMACGGWAIQEj755BN16dJF69evD/QoAEIMsQQgJEycOFG5ubmBHgNACCKWAAAADIglAAAAA2IJQFCbMmWKkpKSlJOTI0kaPny4kpKSNGXKFM+aw4cP6x//+IfuuusuXXvttUpJSVHPnj311ltv6cyZM+fd57hx47Rq1Sp16dJFLVq0UKdOnbRs2TLt379fSUlJ6tChg9xutxYsWKDu3bvr2muv1Y033qgnnnhCWVlZkqSjR49q9OjRuvXWW9WiRQt17NhRo0aNUn5+vn/+gAB4LTzQAwCAN2rXrq3rrrtO27ZtU2FhoRo0aKC4uDjVrl1bkrR582YNGjRIeXl5ioiIUMOGDWXbtrZv365t27Zp6dKlmjNnjq688spz9r1x40a99dZbql69uho3bqysrCw1a9bMs93tdmvw4MFauXKlEhIS1KBBA2VnZ2vVqlXauHGjZs2apcGDB+vQoUNq0KCB6tSpo71792rBggXasWOHFi5cKMuy/PZnBeDiWLZt24EeAgC81alTJ+Xk5GjMmDHq1auXJCk3N1f33nuv8vLy1Lt3bz3//POKiYmRJP3000967rnntHXrVrVu3VoLFizw7GvKlCmaOnWqJKlz586aMGGCKlWqpKNHjyo2Nlb79+/XbbfdJkkKDw/XmDFj1L17d1mWpczMTPXu3VtOp1MOh0NNmzbVxIkT1bBhQ0nSRx99pGHDhkmS5s+fr7Zt2/rrjwjAReJtOAAh680331ReXp46deqk0aNHe0JJkurXr6/p06eratWq2rRpk7766qtS9zF06FBVqlRJkhQbG3vO9vvvv189evTwnCFq0qSJJ6Rs29brr7/uCSVJ6tmzpxITEyVJO3bs8Ml/J4CKRSwBCFmrVq2SJN17772lbq9Zs6bat28vSVqzZs0526+88krVq1fPeIxbb731nOfOxtBVV12l+vXrn7M9Pj5eknTy5EnjvgFcGrhmCUBIOnXqlOei7+nTpys1NbXUdWfXZGdnn7PtbNSYnL026tciIiIklX4m6tfbuQoCCA7EEoCQ9OuzNpmZmRdcX9qn0ypXrnzBn4uKijrvNoeDk/dAKCCWAISkX0fMJ598oiZNmgRwGgDBjF97AISkmJgY1axZU5K0a9eu867LyMhQenq6jh8/7q/RAAQZYglASDj7abRfXwd09uLrd955R263+5yfyc/P16OPPqru3btr3rx5fpkTQPAhlgCEhOjoaEn/uWBbkgYMGKDo6Ght3rxZzz//vI4ePerZlpOTowEDBujYsWOqVq2a+vbt6/eZAQQHrlkCEBKaN2+uzMxMzZkzR19//bU6d+6sQYMGadKkSRoyZIg+/fRTrVixQr/73e9UVFSkPXv2qLi4WNHR0Zo9e7bi4uIC/Z8A4BJFLAEICUOHDpXT6dS6deuUnZ3t+W62jh07atmyZZo7d67Wrl2r3bt3y+VyKTExUe3bt9ef/vSnC95LCcDlja87AQAAMOCaJQAAAANiCQAAwIBYAgAAMCCWAAAADIglAAAAA2IJAADAgFgCAAAwIJYAAAAMiCUAAAADYgkAAMCAWAIAADAglgAAAAyIJQAAAIP/B6fyxwpmSxaOAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Plot\n", "ax = sns.pointplot(data=p,\n", " x='term', y='estimate',\n", " linestyle='none')\n", "# Add errors\n", "ax.vlines(p['term'], p['conf_low'], p['conf_high'])" ] }, { "cell_type": "markdown", "id": "a15aaee3-e020-4cce-957c-fd7d3e9c8616", "metadata": {}, "source": [ "The *p*-value of this difference is significant, as it shows that we'd not expect to see a difference as large as this if the difference between the interventions was zero.\n", "\n", "But that's often a bit of a silly hypothesis to begin with. We'd not run the intervention (or study, or wherever we get our data!) if we thought the difference was really *nothing* to start with. \n", "\n", "Instead, imagine the scenario that your new intervention is out to beat an existing one. That one showed a difference of 7 units by age 2. Is this intervention better than that? Testing this is simple:" ] }, { "cell_type": "code", "execution_count": 8, "id": "9fc0bc6a-da86-40cd-8f58-3d673ff4991e", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "shape: (1, 8)
termestimatestd_errorstatisticp_values_valueconf_lowconf_high
strf64f64f64f64f64f64f64
"b1-b2=6"3.4902941.496982.3315570.0197245.6639050.5562686.42432
" ], "text/plain": [ "shape: (1, 8)\n", "┌─────────┬──────────┬───────────┬──────┬─────────┬──────┬───────┬───────┐\n", "│ Term ┆ Estimate ┆ Std.Error ┆ z ┆ P(>|z|) ┆ S ┆ 2.5% ┆ 97.5% │\n", "│ --- ┆ --- ┆ --- ┆ --- ┆ --- ┆ --- ┆ --- ┆ --- │\n", "│ str ┆ str ┆ str ┆ str ┆ str ┆ str ┆ str ┆ str │\n", "╞═════════╪══════════╪═══════════╪══════╪═════════╪══════╪═══════╪═══════╡\n", "│ b1-b2=6 ┆ 3.49 ┆ 1.5 ┆ 2.33 ┆ 0.0197 ┆ 5.66 ┆ 0.556 ┆ 6.42 │\n", "└─────────┴──────────┴───────────┴──────┴─────────┴──────┴───────┴───────┘\n", "\n", "Columns: term, estimate, std_error, statistic, p_value, s_value, conf_low, conf_high" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Difference\n", "p = me.predictions(cog_mod,\n", " newdata=me.datagrid(cog_mod,\n", " trt=['Y', 'N'],\n", " age=2),\n", " hypothesis='b1-b2 = 7') \n", "# Subtract one from the other, is it equal to 7?\n", "p" ] }, { "cell_type": "markdown", "id": "0c29a935-55af-41d0-ae71-5d5b5c409d7f", "metadata": {}, "source": [ "The predicted difference from seven - *not* zero - is about 2.5, and it is not significant. But note the confidence interval - it might be just less than 7, or quite a bit above it. So while it is not statistically significant, this might well be worth pursuing! This is an example of a large prediction with a wide interval. If the null value is zero this is comfortably significant, but if its 7, it is not - but the interval indicates the effect could be rather a bit larger than 7. Circumstance dictates what the appropriate outcome is.\n", "\n", "You can test almost any hypothesis with `marginaleffects` with a careful application, and often, the hypothesis of zero is not that interesting!" ] }, { "cell_type": "markdown", "id": "55783f2d-ebd3-4e3e-a8aa-68977540735d", "metadata": {}, "source": [ "## Testing an interval hypothesis\n", "We can also express a hypothesis as an interval. \n", "\n", "Imagine this somewhat complex scenario:\n", "- A previous, more expensive intervention produces a change of 7 points.\n", "- We are tasked with seeing if our intervention is *at least as good as* - i.e., *is equivalent* to the existing one, for children aged 2.\n", "- The current intervention will replace the existing one, if it can be shown to be equivalent to give an increase of at least 6.5 - 12.5.\n", "- The developers would like it to be better than 7, but they'll accept anywhere between 6.5 and 12.5 as equivalent.\n", "\n", "Testing this kind of 'interval' hypothesis is easy, but a bit confusing philosophically. A significant result means our estimate is *within* the bounds specified, a non-significant one indicating one or both of the bounds are within the confidence interval and cannot be excluded - i.e., the effect could be lower than 6.5 or bigger than 12.5. We do this using the `equivalence` keyword and a list of the lower and upper bounds, in addition to the standard hypothesis:" ] }, { "cell_type": "code", "execution_count": 9, "id": "b39181dc-a721-480a-98d5-3831f34b4f16", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "shape: (1, 13)
termestimatestd_errorstatisticp_values_valueconf_lowconf_highstatistic_noninfstatistic_nonsupp_value_noninfp_value_nonsupp_value_equiv
strf64f64f64f64f64f64f64f64f64f64f64f64
"b1=b2"9.4902941.496986.3396282.3032e-1032.0156416.55626812.424321.997552-2.0105190.0228830.0221880.022883
" ], "text/plain": [ "shape: (1, 8)\n", "┌───────┬──────────┬───────────┬──────┬─────────┬─────┬──────┬───────┐\n", "│ Term ┆ Estimate ┆ Std.Error ┆ z ┆ P(>|z|) ┆ S ┆ 2.5% ┆ 97.5% │\n", "│ --- ┆ --- ┆ --- ┆ --- ┆ --- ┆ --- ┆ --- ┆ --- │\n", "│ str ┆ str ┆ str ┆ str ┆ str ┆ str ┆ str ┆ str │\n", "╞═══════╪══════════╪═══════════╪══════╪═════════╪═════╪══════╪═══════╡\n", "│ b1=b2 ┆ 9.49 ┆ 1.5 ┆ 6.34 ┆ 2.3e-10 ┆ 32 ┆ 6.56 ┆ 12.4 │\n", "└───────┴──────────┴───────────┴──────┴─────────┴─────┴──────┴───────┘\n", "\n", "Columns: term, estimate, std_error, statistic, p_value, s_value, conf_low, conf_high, statistic_noninf, statistic_nonsup, p_value_noninf, p_value_nonsup, p_value_equiv" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Equivalence\n", "eq = me.predictions(cog_mod,\n", " newdata=me.datagrid(cog_mod,\n", " trt=['Y', 'N'],\n", " age=2),\n", " hypothesis='b1=b2',\n", " equivalence=[6.5, 12.5])\n", "eq" ] }, { "cell_type": "markdown", "id": "a0b7f369-ba20-4765-a580-f1db5eaacf0f", "metadata": {}, "source": [ "The `p_value_equiv` contains the information we need. Illustrating this graphically:" ] }, { "cell_type": "code", "execution_count": 10, "id": "8a0ce2dd-0b63-4179-9245-fd4f92cd22e5", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAksAAAHHCAYAAACvJxw8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAvIklEQVR4nO3de5iN9f7/8de9ZoaZwdCMZjBOsTPIyBRKir4iqV1ItKPDtu3Y6SBXB+WrfB32bpdjzqRiJDr5UmkjUikS+tI4zcSgMdrjOAyzmJm17t8f/azdZHzMWGvWspbn47q6Lmvdn7nvd65r8XSve93Lsm3bFgAAAErkCPQAAAAAlzJiCQAAwIBYAgAAMCCWAAAADIglAAAAA2IJAADAgFgCAAAwCA/0AMEuLS1NhYWFcjgcqlixYqDHAQAApXDmzBm53W5FREQoOTnZuJZY8lJhYaFs25bL5VJ+fn6gxwEAAGVQWFh4wTXEkpccDodcLpcsy1JUVFSgxwEAAKXgdDpl27YcjgtfkUQsealixYrKz89XVFSUmjRpEuhxAABAKezYsUP5+fmluoSGC7wBAAAMiCUAAAADYgkAAMCAWAIAADAglgAAAAyIJQAAAANiCQAAwIBYAgAAMOCmlD7icts67rzwLdMBAEDgudx2qdcSSz6y/ZcT+uO8FYEeAwAAlMJrHePU8IqIUq3lbTgAAAADYgkAAMCAt+F8pGnNGG0Zfn2gxwAAAKWQlfmTzpx2lmotseQjYQ5LVaNK994nAAAIrAMOq9RriSUAKMHpQpemf7nb8/ixWxsqMiIsgBMBCBRiCQBKcKbIrddX/eR5/JebryKWgMsUF3gDAAAYEEsAAAAGxBIAAIABsQQAAGBALAEAABgQSwAAAAbEEgAAgAGxBAAAYEAsAQAAGBBLAAAABsQSAACAAbEEAABgQCwBAAAYEEsAAAAGxBIAAIBBUMbS3r171aJFC40cObLE7U6nU9OmTVPXrl3VokULNW/eXF26dNHYsWOVm5vr32EBAEBQCw/0AGV1+PBhDRw4UE6ns8Ttubm5euihh5SRkaEqVaooJSVFYWFhSktL0xtvvKHPPvtM7777rmrUqOHnyQEAQDAKqjNLO3bsUO/evbV79+7zrhkzZowyMjLUunVrrVixQm+//bZmz56tzz//XLfccouys7P10ksv+XFqAAAQzIIilo4fP64xY8aoV69e2rdvn2rXrl3iutOnT+vTTz+VJP3zn/9UbGysZ1tMTIxeffVVWZalNWvW6NixY36ZHQAABLegiKXU1FTNnj1bsbGxmj59urp161biuiNHjuiaa67Rddddp8TExHO2x8XFqWrVqrJtWwcPHiznqQEAQCgIimuWatSooSFDhqh3796KjIzUtm3bSlyXmJiod99997z72bdvn3Jzc+VwOJSQkFBe4wIAgBASFLHUs2dPn+xn7NixkqSbbrpJ1apV88k+z7JtWy6Xy6f7BBA47t+9nt0ul1yuoDgZD6AUbNsu9dqgiCVfmDFjhlasWKHIyEg9//zzPt+/0+nU5s2bfb5fAIFxqsBd7HFaWpoqVSCWgMvRZRFLkyZN0tSpU+VwOPSPf/xDSUlJgR4JAAAEiZCOpYKCAg0bNkxLlixReHi4XnnlFd11113lcqyoqCgiDAghJ5yF0pJVnsfJycmKiYoI4EQAfCk9Pf2892z8vZCNpSNHjujxxx/X//3f/6ly5cp6/fXXdfPNN5fb8SzLUlhYWLntH4B/OcLcv3scxmscCCGWZZV6bUjG0s8//6w///nPys7OVu3atTVjxgxdffXVgR4LAAAEoZCLpZycHD388MP65Zdf1Lx5c82YMUNxcXGBHgsAAASpkIul5557Tr/88ouSkpI0Z84cVapUKdAjAQCAIBZSsfTtt99q/fr1kqQqVapo+PDh5107aNAg1alTx1+jAQCAIBVSsbR69WrPrzdu3Ghc+8gjjxBLAADggiy7LLewxDl27Nih/Px8RUdHq0mTJoEeB4CPHHcW6toRKzyPtwy/XVW5dQAQMsry9ze3owUAADAglgAAAAyIJQAAAANiCQAAwIBYAgAAMCCWAAAADIglAAAAA2IJAADAgFgCAAAwIJYAAAAMiCUAAAADYgkAAMCAWAIAADAglgAAAAyIJQAAAANiCQAAwIBYAgAAMCCWAAAADIglAAAAA2IJAADAgFgCAAAwIJYAAAAMiCUAAAADYgkAAMCAWAIAADAglgAAAAyIJQAAAANiCQAAwIBYAgAAMCCWAAAADIglAAAAA2IJAADAgFgCAAAwIJYAAAAMiCUAAAADYgkAAMCAWAIAADAglgAAAAyIJQAAAANiCQAAwCCoY2nv3r1q0aKFRo4c6ZN1AAAAvxe0sXT48GENHDhQTqfTJ+sAAABKEpSxtGPHDvXu3Vu7d+/2yToAAIDzCQ/0AGVx/PhxzZo1S6mpqSooKFDt2rW1f//+i14HAABwIUF1Zik1NVWzZ89WbGyspk+frm7dunm1DgAA4EKCKpZq1KihIUOGaPny5erQoYPX6wAAAC4kqN6G69mzp0/X+ZJt23K5XH4/LoDy4f7d69ntcsnlCqp/XwIwsG271GuDKpYuZU6nU5s3bw70GAB85FSBu9jjtLQ0VapALAGXI175AAAABpxZ8pGoqCglJSUFegwAPnLCWSgtWeV5nJycrJioiABOBMCX0tPTS30PRmLJRyzLUlhYWKDHAOAjjjD37x6H8RoHQohlWaVey9twAAAABsQSAACAAbEEAABgQCwBAAAYEEsAAAAGll2WW1jiHDt27FB+fr6io6PVpEmTQI8DwEeOOwt17YgVnsdbht+uqtw6AAgZZfn7mzNLAAAABsQSAACAAbEEAABgQCwBAAAYEEsAAAAGxBIAAIABsQQAAGBALAEAABgQSwAAAAbEEgAAgAGxBAAAYEAsAQAAGBBLAAAABsQSAACAAbEEAABgQCwBAAAYEEsAAAAGxBIAAIABsQQAAGBALAEAABgQSwAAAAbEEgAAgAGxBAAAYEAsAQAAGBBLAAAABsQSAACAAbEEAABgQCwBAAAYEEsAAAAGxBIAAIABsQQAAGBALAEAABgQSwAAAAbhvtrR8ePHtW7dOmVmZiovL09DhgzRmTNntGXLFrVu3dpXhwEAAPArr2PJtm1NnjxZb7/9tk6fPu15fsiQIdq/f78eeeQRtWjRQlOnTlVsbKy3hwMAAPArr9+Ge/755zV9+nQ5nU5Vq1ZNUVFRnm25ubmybVubN2/WQw89JKfT6e3hAAAA/MqrWFqxYoU++eQTxcbG6o033tC6devUuHFjz/brr79e8+fPV1xcnDIzM5Wamur1wAAAAP7kVSy99957sixL48aN0y233FLimuuvv14TJ06Ubdtavny5N4cDAADwO69iaevWrapZs6ZuvPFG47qWLVsqMTFRe/fu9eZwAAAAfudVLOXn56tatWqlWhsbG6uioiJvDgcAAOB3XsVS9erVtW/fPtm2bVxXWFiovXv3qnr16t4cDgAAwO+8iqXWrVsrPz9f7777rnHd3LlzlZeXp5YtW3pzOI+9e/eqRYsWGjly5HnXrF27Vn379lWbNm2UkpKiHj166IMPPrhg2AEAAPyWV7HUt29fORwOvfrqq0pNTdWxY8eKbT9y5IgmTpyo8ePHy+Fw6MEHH/RqWEk6fPiwBg4caLwNwYIFC9S3b19t2LBBTZs21Q033KDdu3dr2LBhGjp0qNczAACAy4dXN6Vs3Lixhg4dqtGjR+uVV17RK6+84tnWpk0b5ebmSvr1xpWDBg1S8+bNvRp2x44dGjRokPbt23feNXv27NGoUaNUuXJlzZs3T02bNpUkHThwQI888ogWLVqk9u3b64477vBqFgAAcHnw+qaUffr00cyZM5WUlCTbtj3/HTt2TLZtq27duho/frwee+yxiz7G8ePHNWbMGPXq1Uv79u1T7dq1z7t29uzZcrlc6tevnyeUJKlWrVp6+eWXPWsAAABKwyffDdeuXTu1a9dO2dnZ+umnn5SXl6eoqChdddVVatiwodf7T01N1ezZs1WjRg0NHz5c27Zt05QpU0pcu3r1aknS7bfffs62m266SVWqVFFaWppycnKUkJDg9WwAACC0eRVLGzZsUJUqVTx37U5MTFRiYmKJa7/55hvt27dPffr0KfNxatSooSFDhqh3796KjIzUtm3bSlx3+PBhHTlyRBEREWrQoME528PCwtSgQQNt2bJF6enpPo0l27blcrl8tj8AgeX+3evZ7XLJ5fL6ZDyAS0RZPvDlVSw99NBDatmypd55550Lrp04ceJFx1LPnj1Lte7gwYOSfr2lgcNR8h9q8fHxxdb6itPp1ObNm326TwCBc6rAXexxWlqaKlUgloDLUalj6eTJk+d82k2STp8+raysrPP+nG3bys7OVmZmZrl/bD8/P1+SFBkZed41FStWLLYWAADApNSxdOrUKd199906c+aM5znLsrRt27YSrw8qSYsWLco8YFmcPZtkWdYF1/o63KKiopSUlOTTfQIInBPOQmnJKs/j5ORkxURFBHAiAL6Unp5uvA3Rb5U6lhISEtS3b19Nnz7d85xlWaWOjlq1amnYsGGlPdxFqVSpkqRfz3adz9nYi46O9umxLctSWFiYT/cJIHAcYe7fPQ7jNQ6EkNKcWDmrTNcsDRw4UPfdd5+kX8/MdOzYUcnJyZo4ceJ5f8bhcCg6OlpVq1Yty6EuytkLtg8fPizbtkv8jTh7rdLZa5cAAABMyhRLERERxT7t1qpVKyUlJZ33E3D+Vq1aNSUkJCgnJ0f79u1T/fr1i213uVzKzMyUJN4yAwAApeLVRzvmzZtX7m+tlVX79u0lSStWrDhn27fffqu8vDw1btxYNWrU8PdoAAAgCPn8c7But7vYfwUFBTpx4oQyMjI0c+ZMXx/uHH369FFYWJhmzpxZ7KP8Bw4c0KhRoyRJAwYMKPc5AABAaPD6Dt5ff/21Jk+erIyMDBUUFFxwfXmHSuPGjfX0009r3Lhx6t27t1q3bq3IyEitX79e+fn56tmzp+68885ynQEAAIQOr2Jp69atGjhwoFwu1wU/FRceHq6UlBRvDldq/fv3V8OGDTVnzhz9+OOPsixLDRs21AMPPKDu3bv7ZQYAABAavIql1NRUFRUV6eqrr9Zf//pXRUZGatCgQercubPuv/9+/fvf/9aiRYu0ceNGtWzZUnPmzPHJ0E8++aSefPJJ45rbbrtNt912m0+OBwAALl9exdLGjRsVFhamSZMm6aqrrpIk1axZU1lZWbrpppskSd27d9dTTz2llStXaunSpbrrrru8nxoAAMBPvLrA+8iRI6pVq5YnlKRfrxn67fVLlmXpxRdflCQtWrTIm8MBAAD4ndefhqtWrVqxx/Xr15fL5dKePXs8z9WqVUv16tVTenq6t4cDAADwK69iKS4uTocOHSr2XO3atSVJu3btKvZ8pUqVlJub683hAAAA/M6rWGrWrJlycnK0Zs0az3MNGjSQbdvasGGD57nTp0/r559/9stXngAAAPiSV7HUtWtX2batp556SmPGjFFRUZFSUlJUtWpVffjhh1q8eLEyMjI0bNgw5eXlFbu2CQAAIBh4FUu33Xab7rrrLjmdTs2dO1dhYWGKjIzUAw88oKKiIr344ovq2rWrli5dKsuy9PDDD/tqbgAAAL/w+g7e48aNU9u2bfXtt9/KsixJ0hNPPKGcnBwtXrxYtm0rLCxMffv21e233+71wAAAAP7kdSxJ0r333qt77733PzsND9crr7yiwYMH68CBA6pbt65iY2N9cSgAAAC/8kksnU98fLzi4+PL8xAAAADlyiexVFBQoMzMTJ08efKC3xHXqlUrXxwSAMrNwROn9f7GrGLPzVu3V71a1lF8TGSApgIQKF7H0oQJE5SamqrTp09fcK1lWdq+fbu3hwSAcpG2/7hmfL1by7f+W0Xu4v/wG7siQxNX/qTOzWrob+0aKrk2t0IBLhdexdKcOXM0c+bMUq+/0FknAAiUd77bp5eXbJXb8MdUkdvW0h9/0b/SftHIrs304I31/DcggIDxKpbef/99WZale++9V48//rji4+MVHl6ul0EBgM+9890+DVu8tdTr3bY0bPFWWZbU5waCCQh1Xt1nKSsrS3FxcRo1apRq1apFKAEIOmn7j+vlJaUPpd96afFWpe0/7uOJAFxqvIqlypUr68orr5TD4fX38QJAQMz4erfxrTcTty3N/Hq3bwcCcMnxqnJatWqlPXv26NSpU76aBwD85uCJ01q+9d9e7WPZ1n/r4IkLf8AFQPDyKpYGDhwol8ulESNGcPE2gKDz8ZYD53zqrayK3LY+3nLARxMBuBR5dZFR48aN9frrr+uJJ57Q5s2b1bZtW8XGxnq+9qQkTzzxhDeHBACf2X3opI/2w9l1IJR5FUsnT57UrFmz5Ha7lZWVpYULF17wZ4glAJeKU2dcPtpPkU/2A+DS5FUsTZgwQZs3b5YkXXnllapVq5YiIiJ8MRcAlLtKFcN8tB8+CQyEMq9e4atWrZJlWfrHP/6h7t27+2omAPCLhldW9tF+KvlkPwAuTV5d4H306FHVq1ePUAIQlO65tpbCHee/xrI0wh2W7rm2lo8mAnAp8iqW4uPjFRbmm9PYAOBv8TGR6tyshlf7uKNZDb5cFwhxXsVS586dlZmZqW3btvlqHgDwq7+1a6iLPbnksKQB7Rr6diAAlxyvYumxxx5T3bp19dhjj2nZsmXcnBJA0EmuXVUjuza7qJ8d1a2ZkmtX9fFEAC41Xl3gPXLkSNWpU0fffPONBg8eLMuyVKVKFUVFRZW43rIsrV692ptDAoDPPXhjPVnWr9/1Vpp7VDqsX0OJL9EFLg9exdLHH3/s+bVt27JtW8ePH9fx4yV/saTpZpUAEEh9bqin5onVNPPr3Vq29d8l3tk73GHpjmY1NKBdQ84oAZcRr2LplVde8dUcABBwybWrakrv63TwxGm9vzFLY1dkeLY9e3sj9WpVR/FVuJgbuNx4FUvcMgBAKIqPidRDbeoXi6WH2tRX1Shuugtcjry6wBsAACDUlfrM0rp16yRJ1113nSpWrFjsubJo06ZNmX8GAAAgUEodS3379pXD4dDSpUt11VVXeZ4ry0XblmVp+/btZZ8SAAAgQMp0zZLb7T7nOdsuxedsL2ItAADApaDUsbRz585SPQcAABBK/HqB95EjR/x5OAAAAK95FUu33XabBg8eXKq1f/rTn7jVAAAACDpexVJ2drYOHjx4wXVut1uHDh3SsWPHvDkcAACA35X6mqVdu3Zp+PDh5zyfkZGhPn36nPfnbNtWTk6ODhw4oFq1al3clAAAAAFS6lj6wx/+oMjISH377bee5yzLUl5enjZt2lSqfTz44INlnxAAACCAynTrgJdeekmffvqp5/GUKVNUq1Yt3Xvvvef9GcuyVKlSJTVp0kQ33HDDxU9aRv/617/07rvvavv27bJtW/Xr11ePHj3Us2dPVahQwW9zAACA4FamWKpfv76eeOIJz+MpU6aoZs2axZ67FAwfPlwLFy6UJDVq1Ei1a9fWzp07NXLkSH322WeaNm2aqlblG8MBAMCFefVFuqtWrfJ89cmlYsmSJVq4cKEiIiI0ZswYdenSRZJUWFioV199VfPmzdPo0aM1ZsyYAE8KAACCgVefhktMTFT16tVL3Hb69Gl98cUXWrlypXJzc705TJksWLBAkvToo496QkmSIiIi9MILL6hhw4b6+OOPlZGRcb5dAAAAeHh1ZkmScnJyNH36dNWqVUv9+/eXJO3evVt9+/bVoUOHJElRUVEaPXq07rzzTm8Pd0Hp6emSpI4dO56zLTw8XK1atdLu3bv11VdfqVGjRuU+DwAACG5exdLRo0fVq1cvHTx4ULfeeqvn+ZdfflkHDx70XNx98uRJPf/880pKSlLDhg29ndnI5XJJkipXrlzi9vDwX/+Xd+/e7dPj2rbtOTaA4Of+3evZ7XLJ5fLrlx4AKEdl+b5ar2Jp7ty5ysnJUb169XT//fdLkvbt26dNmzYpLCxM8+fPV4sWLTR+/HjNmjVLc+bM0ahRo7w55AU1aNBAO3bs0Pfff6969eoV22bbtn744QdJvv/qFafTqc2bN/t0nwAC51RB8S8OT0tLU6UKxBJwOfLqlf/1118rPDxcb775pufM0pdffilJuu6669SiRQtJ0pNPPqmYmBh999133hyuVM7exmDs2LHF4sXtdmvSpEnavn27JKmgoKDcZwEAAMHPqzNLWVlZql+/vmrXru15bu3atbIsSzfddJPnuYiICNWuXdvnb32VpE+fPvruu++0atUqPfDAA0pOTlb16tWVnp6unJwc/elPf9LChQs9b8f5SlRUlJKSkny6TwCBc8JZKC1Z5XmcnJysmKiIAE4EwJfS09PldDpLtdarYnC5XMVu8FhUVKQNGzZIklq3bl1srdPplGVZ3hyuVMLCwjRlyhS9++67+uCDD7Rjxw5FR0frhhtu0JQpU7Rnzx4tXLjQ5/dZsixLYWFhPt0ngMBxhLl/9ziM1zgQQsrSJF7FUmJiorKzs1VYWKiIiAht2LBB+fn5qly5suctOOnXT8xlZWWpTp063hyu1BwOhx588MESv17l888/98wOAABwIV5ds5ScnKwTJ05o7Nix2rlzpyZOnCjLstS+fXvPv8COHDmi5557Ti6XS23atPHJ0CY///yzvvnmG89tC35v7dq1kqTmzZuX+ywAACD4eRVLjz76qCIjI5Wamqru3btry5YtCgsL06OPPipJ2rhxo9q3b68NGzaoSpUq+stf/uKToU0++ugj9evXTx9++OE527Zv367NmzerWrVqatu2bbnPAgAAgp9XsdSgQQO99dZbSk5OVoUKFdSoUSNNnz5djRs3liTFx8erqKhIV199tRYsWFDsQvDy0rFjR1mWpTlz5igrK8vz/C+//KJnn31Wtm1rwIABio6OLvdZAABA8PP6I2EpKSl6//33S9xWu3ZtLV682BNP/pCcnKy//vWveuONN3T33XerVatWkqT169frzJkz6t69u/785z/7bR4AABDcfPb5+TVr1mj16tXKzMxUXl6ePvroI508eVIrV65UQkKCrrjiCl8d6oKeeeYZ1alTRwsWLNB3332nSpUqqXnz5urdu7e6dOnil0/lAQCA0OB1LB05ckRPP/20Nm7cKOnXu2SfjZEDBw5oypQpmjdvnmbNmqVrr73W28OVimVZuv/++z13FQcAALhYXl2zVFBQoH79+mnDhg2qVKmSOnXqpISEhP/s3OFQtWrVdPz4cfXt21fZ2dleDwwAAOBPXsXS/PnztXPnTrVo0UIrVqzQpEmTit2/qFGjRlq5cqVSUlLkdDr19ttvez0wAACAP3kVS0uXLpXD4dCYMWMUGxtb4prKlStr7NixCgsL05o1a7w5HAAAgN95FUuZmZlq2LDhBe/MnZiYqPr16+uXX37x5nAAAAB+51Usud3uCy/6/yIiIvheJQAAEHS8iqXExETt3btXJ0+eNK47duyYfvrpJ76PDQAABB2vYql9+/YqLCzUmDFjjOtGjx4tl8ulW265xZvDAQAA+J1X91nq16+fPvroI73//vs6cuSI7r77buXl5UmSdu/erYyMDM2fP1+bNm1SpUqVuHM2AAAIOl7FUlxcnKZNm6aBAwdq5cqVWrVqlWfbH//4R0m/3qQyOjpa48ePL3YPJgAAgGDg1dtwknT99dfr448/1sMPP6yaNWvKtm3Pf3Fxcbrvvvu0ePFitWvXzhfzAgAA+JVPvhsuISFBQ4cO1dChQ5Wfn6+8vDxFR0erSpUqvtg9AABAwPjsi3TPio6OVnR0tK93CwAAEBBevw0HAAAQyoglAAAAA2IJAADAgFgCAAAwIJYAAAAMiCUAAAADYgkAAMCAWAIAADAglgAAAAyIJQAAAANiCQAAwIBYAgAAMCCWAAAADIglAAAAA2IJAADAgFgCAAAwIJYAAAAMiCUAAAADYgkAAMCAWAIAADAglgAAAAyIJQAAAANiCQAAwIBYAgAAMCCWAAAADIglAAAAA2IJAADAgFgCAAAwCA/0AOXpq6++0pw5c5SWlqbTp08rISFBt956qx577DFVr1490OMBAIAgELJnlt5++231799f69at09VXX6327duroKBA77zzjrp166aff/450CMCAIAgEJKxlJ2drXHjxqlChQqaO3euFixYoKlTp2rVqlXq0qWLDh06pL///e+BHhMAAASBkIyldevWqbCwUG3bttUNN9zgeb5ChQp6+umnJUnr168P0HQAACCYhGQshYWFSZIOHjx4zrbDhw9Lkq644gq/zgQAAIJTSMZSmzZtFBERoW3btmn48OHKzs6W0+nUunXr9OKLL0qS+vfvH+ApAQBAMAjJT8PVqFFDEyZM0H//939r4cKFWrhwoWfbFVdcoSlTpqhTp04+PaZt23K5XD7dJ4DAcf/u9ex2ueRyheS/L4HLkm3bpV4bkrEkSU2bNtXtt9+u//3f/1VycrKqVaumbdu26eDBg5o1a5YaN26sOnXq+Ox4TqdTmzdv9tn+AATWqQJ3scdpaWmqVIFYAi5HIRlLO3bsUN++fVWxYkV98MEHatq0qSSpsLBQ48eP11tvvaWHHnpIn332maKjowM8LQAAuJSFZCyNHj1ax44d09SpUz2hJEkRERF6/vnntWXLFm3atEkffvihHn74YZ8cMyoqSklJST7ZF4DAO+EslJas8jxOTk5WTFREACcC4Evp6elyOp2lWhtysXTmzBn98MMPsixLbdu2PWe7ZVlq3769Nm3apK1bt/rsuJZleT6FByD4OcLcv3scxmscCCGWZZV6bci9AX/ixAm53W5jvJx9vqioyJ+jAQCAIBRysRQXF6dq1arJ7Xbryy+/LHHNt99+K0lq0qSJHycDAADBKORiyeFw6IEHHpAk/f3vf1dGRoZnm9vt1pQpU7R27VrFxMSoR48egRoTAAAEiZC7ZkmSHn/8ce3cuVOrV69W165ddd1116lq1arauXOnsrOzFR0drddff12xsbGBHhUAAFziQjKWIiIiNH36dC1atEiLFi3Szp07debMGcXHx6tXr1569NFHVbdu3UCPCQAAgkBIxpL061XuPXr04K02AADglZC7ZgkAAMCXiCUAAAADYgkAAMCAWAIAADAglgAAAAyIJQAAAANiCQAAwIBYAgAAMCCWAAAADIglAAAAA2IJAADAgFgCAAAwIJYAAAAMiCUAAAADYgkAAMCAWAIAADAglgAAAAyIJQAAAANiCQAAwIBYAgAAMCCWAAAADIglAAAAA2IJAADAgFgCAAAwIJYAAAAMiCUAAAADYgkAAMCAWAIAADAglgAAAAyIJQAAAANiCQAAwIBYAgAAMCCWAAAADIglAAAAA2IJAADAgFgCAAAwIJYAAAAMiCUAAAADYgkAAMCAWAIAADAID/QAvtahQwdlZ2dfcF3r1q01b948P0wEAACCWcjFUseOHXX06NESt9m2rWXLlqmoqEjXXHONnycDAADBKORiaejQoefdNnXqVBUVFalVq1Z69tln/TgVAAAIVpfNNUvr16/XlClTFBMTo3Hjxik8POQ6EQAAlIPLIpYKCgr08ssvy+126/nnn1dCQkKgRwIAAEHisoilOXPmaO/evUpOTtZ9990X6HEAAEAQCfn3ok6ePKlZs2ZJkp566ilZllUux7FtWy6Xq1z2DcD/3L97PbtdLrlcl8W/L4HLgm3bpV4b8rG0cOFC5eXl6ZprrlG7du3K7ThOp1ObN28ut/0D8K9TBe5ij9PS0lSpArEEXI5C+pXvcrmUmpoqSerfv3+ApwEAAMEopM8sff/998rJyVHVqlXVoUOHcj1WVFSUkpKSyvUYAPznhLNQWrLK8zg5OVkxUREBnAiAL6Wnp8vpdJZqbUjH0vLlyyVJnTt3VoUKFcr1WJZlKSwsrFyPAcB/HGHu3z0O4zUOhJCyXMMc0m/DffXVV5KkO+64I8CTAACAYBWysXTo0CEdOHBA4eHhSklJCfQ4AAAgSIVsLP3444+SpEaNGik6OjrA0wAAgGAVsrGUlZUlSapTp06AJwEAAMEsZGPp2LFjkqSaNWsGeBIAABDMQvbTcIMHD9bgwYMDPQYAAAhyIXtmCQAAwBeIJQAAAANiCQAAwIBYAgAAMCCWAAAADIglAAAAA2IJAADAgFgCAAAwIJYAAAAMiCUAAAADYgkAAMCAWAIAADAglgAAAAyIJQAAAANiCQAAwIBYAgAAMAgP9AAAcCmqGO7QoNuuLvYYwOWJWAKAEkRGhGlwp0aBHgPAJYBY8hGX29ZxZ2GgxwAAAKXgctulXkss+cj2X07oj/NWBHoMAABQCq91jFPDKyJKtZY34QEAAAyIJQAAAAPehvORpjVjtGX49YEeAwAAlEJW5k86c9pZqrXEko+EOSxVjSrde58AACCwDjisUq/lbTgAAAADYgkAAMCAWAIAADAglgAAAAyIJQAAAANiCQAAwIBYAgAAMCCWAAAADLgppZfOnDkjSXI6ndqxY0eApwEAAKXhdP569+6zf4+bEEtecrvdkiTbtpWfnx/gaQAAQFmc/XvchFjyUkREhAoLC+VwOFSxYsVAjwMAAErhzJkzcrvdioi48FeVWbZt236YCQAAIChxgTcAAIABsQQAAGBALAEAABgQSwAAAAbEEgAAgAGxBAAAYEAsAQAAGBBLAAAABsQSAACAAbEEAABgQCwBAAAYEEsAAAAGxBIAAIABsQQAAGBALAEAABgQSwCCwuTJk5WUlKSRI0de9D6WLFmipKQkLVu2zGdzvfDCC0pKStKbb75Zpp87evSoXnvtNd15551q3ry5WrRooa5du2rGjBk6ffq0z+YD4L3wQA8AAP7www8/aMSIEYEeQ5KUlZWlPn36KCcnR7GxsbrxxhtVUFCgLVu2aMKECVqxYoVSU1NVuXLlQI8KQJxZAnAZ+PTTT9WvXz+dOnUq0KNIkl566SXl5OSoS5cuWrVqlWbNmqU5c+Zo2bJluuaaa7Rt2zaNHTs20GMC+P+IJQAha8+ePXrqqaf0zDPPSJKqV68e4Imk/fv3a926dapSpYpGjx6t6Ohoz7aEhAT9z//8j6RfA8+27QBNCeC3eBsOQNBZu3atJk+erO3btysyMlKtWrXSgAEDlJycXGzdyy+/rO+//14pKSkaPXq0RowYocOHD5e4z/Xr1+vhhx8u1fG7d++uf/7zn+c8v3TpUs2ePVu7du1STEyM2rZtq4EDB6p+/fqeNUeOHFFKSooSEhJKfJutQYMGkqS8vDydOnWKt+KASwCxBCCorF27VgsWLFCNGjXUvn17ZWVl6fPPP9fq1as1ceJEderUybO2WbNmevDBB3X77bfLsizjfqtXr6677767VDOkpKSc89yiRYu0a9cuXXXVVfqv//ovZWRkaMmSJfr888/11ltveX7m2muv1cKFC8+77x9//FGSFBMTQygBlwhiCUBQ2bNnj7p3766RI0eqQoUKkqT58+dr5MiRGjp0qFq2bKkrrrhCkjRkyJBS77dhw4ZeXSe0a9cu/e1vf9PTTz8ty7Jk27bGjRunN954Q88++6yWLVumiIgI4z5cLpcmTJggSerSpctFzwLAt7hmCUBQiYuL00svveQJJUnq06eP2rdvrxMnTmjJkiUBmatRo0YaNGiQ5wyWZVl65pln1KhRI+3fv19ffvml8edt29aIESP0448/KjY2Vk888YQfpgZQGsQSgKDSqVMnVapU6ZznO3bsKEn67rvv/D2SJOnuu++Ww1H8j1TLstShQwdJv14TdT5FRUUaOnSo3nvvPVWsWFGTJk1SfHx8uc4LoPR4Gw5AUKlTp06Jz9esWVOSdPDgwYvar7cXeF9orpycnBK3nzhxQoMHD9Y333yj6OhoTZs2Ta1atSrD5ADKG7EEIKhUrFixxOfPfsz+QtcFnY+3F3hfzFz79u3TgAEDtGfPHsXHx2vGjBm65ppryjA1AH8glgAElfOdodm/f7+k/5zJKStvL/Au61w//vijHn30UeXm5qpp06aaPn26atSocdHHB1B+uGYJQFBZs2ZNiTdr/Ne//iVJuvHGG/09kiTp66+/Pue5oqIirVy5UlLxudLT09WvXz/l5ubq1ltv1fz58wkl4BJGLAEIKjt37tTEiROLBdPMmTP1/fffKyEhQffcc09A5vriiy+K3T+pqKhIo0eP1t69e9W0aVPdfPPNkqSCggI9/fTTOnHihG6++WZNnTq12F28AVx6eBsOQFBJSUnRjBkztHz5ciUlJWnXrl3atWuXqlSpokmTJgUsPFJSUjR8+HC99957qlu3rrZu3ar9+/crISFBEyZM8NxSYPHixcrMzJQkORwOvfDCC+fd54gRI0r85B8A/yKWAASVP/7xj+rbt69mzpypL774QpUrV1bXrl315JNPnvcTaf4wYMAA5eTkaO7cuVq1apViY2PVu3dvPf7448W+k2716tWeX5f01t1vDR06lFgCLgGWzTc1AgAAnBfXLAEAABgQSwAAAAbEEgAAgAGxBAAAYEAsAQAAGBBLAAAABsQSAACAAbEEAABgQCwBAAAYEEsAQobT6dT+/fsDPQaAEEMsAQgJn3zyiTp37qx169YFehQAIYZYAhASJkyYoJycnECPASAEEUsAAAAGxBIAAIABsQQgqE2ePFlJSUnKzs6WJA0bNkxJSUmaPHmyZ83hw4f12muv6c4779S1116rlJQU9ejRQ2+99ZbOnDlz3n2OHTtWK1euVOfOndWsWTN16NBBS5cu1f79+5WUlKR27drJ7XZr/vz56tatm6699lrdeOONevzxx7V7925J0tGjRzVq1Cjdeuutatasmdq3b6+RI0cqLy/PP79BALwWHugBAMAbNWvW1HXXXaetW7eqoKBA9erVU1xcnGrWrClJ2rRpkwYOHKjc3FxFRESofv36sm1b27Zt09atW7VkyRLNnj1bV1555Tn73rBhg9566y1VrVpVDRs21O7du9WkSRPPdrfbrUGDBmnFihVKSEhQvXr1lJmZqZUrV2rDhg2aOXOmBg0apEOHDqlevXqqVauW9u3bp/nz52v79u1asGCBLMvy2+8VgItj2bZtB3oIAPBWhw4dlJ2drdGjR6tnz56SpJycHN1zzz3Kzc1Vr1699NxzzykmJkaS9PPPP+vZZ5/Vli1b1LJlS82fP9+zr8mTJ2vKlCmSpE6dOmn8+PGqUKGCjh49qtjYWO3fv1+33XabJCk8PFyjR49Wt27dZFmWMjIy1KtXLzmdTjkcDjVu3FgTJkxQ/fr1JUkfffSRhg4dKkmaN2+eWrdu7a/fIgAXibfhAISsN998U7m5uerQoYNGjRrlCSVJqlu3rqZNm6bKlStr48aN+uqrr0rcx5AhQ1ShQgVJUmxs7Dnb77vvPnXv3t1zhqhRo0aekLJtW6+//ronlCSpR48eSkxMlCRt377dJ/+fAMoXsQQgZK1cuVKSdM8995S4vXr16mrbtq0kafXq1edsv/LKK1WnTh3jMW699dZznjsbQ1dddZXq1q17zvb4+HhJ0smTJ437BnBp4JolACHp1KlTnou+p02bptTU1BLXnV2TmZl5zrazUWNy9tqo34qIiJBU8pmo327nKgggOBBLAELSb8/aZGRkXHB9SZ9Oq1ix4gV/Lioq6rzbHA5O3gOhgFgCEJJ+GzGffPKJGjVqFMBpAAQz/tkDICTFxMSoevXqkqRdu3add116erp27Nih48eP+2s0AEGGWAIQEs5+Gu231wGdvfj6nXfekdvtPudn8vLy9Mgjj6hbt26aO3euX+YEEHyIJQAhITo6WtJ/LtiWpP79+ys6OlqbNm3Sc889p6NHj3q2ZWdnq3///jp27JiqVKmiPn36+H1mAMGBa5YAhISmTZsqIyNDs2fP1tdff61OnTpp4MCBmjhxogYPHqxPP/1Uy5cv1x/+8AcVFhZq7969KioqUnR0tGbNmqW4uLhA/y8AuEQRSwBCwpAhQ+R0OrV27VplZmZ6vputffv2Wrp0qebMmaM1a9Zoz549crlcSkxMVNu2bfWXv/zlgvdSAnB54+tOAAAADLhmCQAAwIBYAgAAMCCWAAAADIglAAAAA2IJAADAgFgCAAAwIJYAAAAMiCUAAAADYgkAAMCAWAIAADAglgAAAAyIJQAAAANiCQAAwOD/Abve9HQzbyAlAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Plot\n", "ax = sns.pointplot(data=eq, x='term', y='estimate')\n", "ax.vlines(eq['term'], eq['conf_low'], eq['conf_high'])\n", "ax.axhline(6.5)\n", "ax.axhline(12.5)" ] }, { "cell_type": "markdown", "id": "9ec3c04c-36d6-4fd3-9f04-4661699161ac", "metadata": {}, "source": [ "This technique can be very useful. Sometimes you may see a small, precise, and significant result. Rather than taking this as evidence, we may test whether the result is within an interval of effects we could consider either important or unimportant (the definition is up to you!), and test *that*. \n", "\n", "The downside is that defining these intervals is difficult, requires justification, and can be frowned on by less statistically literate folk." ] }, { "cell_type": "markdown", "id": "1dce75fe-fd73-4dbf-85a4-79c6ff033b06", "metadata": {}, "source": [ "## One more example\n", "Lets consider another example, revisiting the beauty and wages dataset. First, we read it in and fit a simple model where we predict hourly wages from a `looks` variable that goes from 1 to 5, and an `exper` variable that counts the persons years of experience in the workforce. This shows a statistically significant effect of looks and experience:" ] }, { "cell_type": "code", "execution_count": 11, "id": "9b171329-f0f5-4f0b-a68d-30ef8043df53", "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
OLS Regression Results
Dep. Variable: wage R-squared: 0.064
Model: OLS Adj. R-squared: 0.062
No. Observations: 1260 F-statistic: 42.70
Covariance Type: nonrobust Prob (F-statistic): 1.15e-18
\n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "\n", " \n", "\n", "
coef std err t P>|t| [0.025 0.975]
Intercept 2.5094 0.671 3.742 0.000 1.194 3.825
exper 0.0971 0.011 9.018 0.000 0.076 0.118
looks 0.6373 0.188 3.390 0.001 0.268 1.006


Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified." ], "text/latex": [ "\\begin{center}\n", "\\begin{tabular}{lclc}\n", "\\toprule\n", "\\textbf{Dep. Variable:} & wage & \\textbf{ R-squared: } & 0.064 \\\\\n", "\\textbf{Model:} & OLS & \\textbf{ Adj. R-squared: } & 0.062 \\\\\n", "\\textbf{No. Observations:} & 1260 & \\textbf{ F-statistic: } & 42.70 \\\\\n", "\\textbf{Covariance Type:} & nonrobust & \\textbf{ Prob (F-statistic):} & 1.15e-18 \\\\\n", "\\bottomrule\n", "\\end{tabular}\n", "\\begin{tabular}{lcccccc}\n", " & \\textbf{coef} & \\textbf{std err} & \\textbf{t} & \\textbf{P$> |$t$|$} & \\textbf{[0.025} & \\textbf{0.975]} \\\\\n", "\\midrule\n", "\\textbf{Intercept} & 2.5094 & 0.671 & 3.742 & 0.000 & 1.194 & 3.825 \\\\\n", "\\textbf{exper} & 0.0971 & 0.011 & 9.018 & 0.000 & 0.076 & 0.118 \\\\\n", "\\textbf{looks} & 0.6373 & 0.188 & 3.390 & 0.001 & 0.268 & 1.006 \\\\\n", "\\bottomrule\n", "\\end{tabular}\n", "%\\caption{OLS Regression Results}\n", "\\end{center}\n", "\n", "Notes: \\newline\n", " [1] Standard Errors assume that the covariance matrix of the errors is correctly specified." ], "text/plain": [ "\n", "\"\"\"\n", " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: wage R-squared: 0.064\n", "Model: OLS Adj. R-squared: 0.062\n", "No. Observations: 1260 F-statistic: 42.70\n", "Covariance Type: nonrobust Prob (F-statistic): 1.15e-18\n", "==============================================================================\n", " coef std err t P>|t| [0.025 0.975]\n", "------------------------------------------------------------------------------\n", "Intercept 2.5094 0.671 3.742 0.000 1.194 3.825\n", "exper 0.0971 0.011 9.018 0.000 0.076 0.118\n", "looks 0.6373 0.188 3.390 0.001 0.268 1.006\n", "==============================================================================\n", "\n", "Notes:\n", "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n", "\"\"\"" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Read in looks\n", "looks = pd.read_csv('https://vincentarelbundock.github.io/Rdatasets/csv/wooldridge/beauty.csv')\n", "\n", "# Fit a simple model and examine output\n", "looks_mod = smf.ols('wage ~ exper + looks', data=looks).fit()\n", "looks_mod.summary(slim=True)" ] }, { "cell_type": "markdown", "id": "0fe3b330-9219-4b40-80d8-c6b211dd5d7d", "metadata": {}, "source": [ "As a silly example, lets assume taking up regular exercise, clean eating, and expensive beauty products will shift your looks from a 3 to a 4. Averaging over your years of experience, how much does going from a 3 to a 4 increase your hourly wage?" ] }, { "cell_type": "code", "execution_count": 12, "id": "a3fcb25a-4236-4590-883b-350d3590b360", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "shape: (2, 8)
looksestimatestd_errorstatisticp_values_valueconf_lowconf_high
i64f64f64f64f64f64f64f64
36.3624360.13248148.0253070.0inf6.1027786.622094
46.9997050.20222434.613560.0inf6.6033527.396057
" ], "text/plain": [ "shape: (2, 8)\n", "┌───────┬──────────┬───────────┬──────┬─────────┬─────┬──────┬───────┐\n", "│ looks ┆ Estimate ┆ Std.Error ┆ z ┆ P(>|z|) ┆ S ┆ 2.5% ┆ 97.5% │\n", "│ --- ┆ --- ┆ --- ┆ --- ┆ --- ┆ --- ┆ --- ┆ --- │\n", "│ str ┆ str ┆ str ┆ str ┆ str ┆ str ┆ str ┆ str │\n", "╞═══════╪══════════╪═══════════╪══════╪═════════╪═════╪══════╪═══════╡\n", "│ 3 ┆ 6.36 ┆ 0.132 ┆ 48 ┆ 0 ┆ inf ┆ 6.1 ┆ 6.62 │\n", "│ 4 ┆ 7 ┆ 0.202 ┆ 34.6 ┆ 0 ┆ inf ┆ 6.6 ┆ 7.4 │\n", "└───────┴──────────┴───────────┴──────┴─────────┴─────┴──────┴───────┘\n", "\n", "Columns: looks, estimate, std_error, statistic, p_value, s_value, conf_low, conf_high" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Make a datagrid\n", "dg = me.datagrid(looks_mod,\n", " looks=[3, 4],\n", " exper=np.arange(1, 40)\n", " )\n", "\n", "# Make the predictions\n", "p = me.predictions(looks_mod,\n", " newdata=dg,\n", " by='looks') # average over experience\n", "\n", "p" ] }, { "cell_type": "markdown", "id": "e819cbbe-4a2d-4094-8201-63976aac5242", "metadata": {}, "source": [ "With looks of about 3, you earn 6.36 an hour. With a four, you earn basically 7. Is this difference statistically significant?" ] }, { "cell_type": "code", "execution_count": 13, "id": "423dda04-2807-4e46-b085-ac09745d6fdd", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "shape: (1, 8)
termestimatestd_errorstatisticp_values_valueconf_lowconf_high
strf64f64f64f64f64f64f64
"b2-b1=0"0.6372690.1880083.3895850.000710.4803860.268781.005757
" ], "text/plain": [ "shape: (1, 8)\n", "┌─────────┬──────────┬───────────┬──────┬─────────┬──────┬───────┬───────┐\n", "│ Term ┆ Estimate ┆ Std.Error ┆ z ┆ P(>|z|) ┆ S ┆ 2.5% ┆ 97.5% │\n", "│ --- ┆ --- ┆ --- ┆ --- ┆ --- ┆ --- ┆ --- ┆ --- │\n", "│ str ┆ str ┆ str ┆ str ┆ str ┆ str ┆ str ┆ str │\n", "╞═════════╪══════════╪═══════════╪══════╪═════════╪══════╪═══════╪═══════╡\n", "│ b2-b1=0 ┆ 0.637 ┆ 0.188 ┆ 3.39 ┆ 0.0007 ┆ 10.5 ┆ 0.269 ┆ 1.01 │\n", "└─────────┴──────────┴───────────┴──────┴─────────┴──────┴───────┴───────┘\n", "\n", "Columns: term, estimate, std_error, statistic, p_value, s_value, conf_low, conf_high" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Include hypothesis of zero difference, test to reject\n", "me.predictions(looks_mod,\n", " newdata=dg,\n", " by='looks',\n", " hypothesis='b2 - b1 = 0') # Same as b2=b1 or 4 - 3" ] }, { "cell_type": "markdown", "id": "6d89caee-106b-4c8a-b9fc-42ea873569d6", "metadata": {}, "source": [ "The difference of about 0.64 is significant. Hold on though, all that exercise, expensive products and good food is a LOT of effort. Will it earn you at least an extra 0.50 an hour?" ] }, { "cell_type": "code", "execution_count": 14, "id": "c3229f4c-4041-4af6-8eaf-8103dec5b53a", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "shape: (1, 8)
termestimatestd_errorstatisticp_values_valueconf_lowconf_high
strf64f64f64f64f64f64f64
"b2-b1=0.50"0.1372690.1880080.7301220.4653161.103719-0.231220.505757
" ], "text/plain": [ "shape: (1, 8)\n", "┌────────────┬──────────┬───────────┬──────┬─────────┬─────┬────────┬───────┐\n", "│ Term ┆ Estimate ┆ Std.Error ┆ z ┆ P(>|z|) ┆ S ┆ 2.5% ┆ 97.5% │\n", "│ --- ┆ --- ┆ --- ┆ --- ┆ --- ┆ --- ┆ --- ┆ --- │\n", "│ str ┆ str ┆ str ┆ str ┆ str ┆ str ┆ str ┆ str │\n", "╞════════════╪══════════╪═══════════╪══════╪═════════╪═════╪════════╪═══════╡\n", "│ b2-b1=0.50 ┆ 0.137 ┆ 0.188 ┆ 0.73 ┆ 0.465 ┆ 1.1 ┆ -0.231 ┆ 0.506 │\n", "└────────────┴──────────┴───────────┴──────┴─────────┴─────┴────────┴───────┘\n", "\n", "Columns: term, estimate, std_error, statistic, p_value, s_value, conf_low, conf_high" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Include hypothesis of a 50-cent difference, test to reject\n", "me.predictions(looks_mod,\n", " newdata=dg,\n", " by='looks',\n", " hypothesis='b2 - b1 = 0.50') " ] }, { "cell_type": "markdown", "id": "8f9daf70-378a-44aa-bb94-83eed9ae234b", "metadata": {}, "source": [ "This difference is not statistically significant, so you conclude the conclude there's no evidence it will improve your wages above 50 cents an hour.\n", "\n", "Imagine now you say - this increase of a 3 to a 4 is going to be a lot of work.\n", "\n", "Anything from an increase of 0.01 to 1 dollar is absolutely not worth it, and so I'll consider that increase equivalent to what I currently earn - all the gym memberships and beauty products will cost me more than that. An increase of 0.01 to 1 dollar is my interval hypothesis.\n", "\n", "You can test that hypothesis like so:" ] }, { "cell_type": "code", "execution_count": 15, "id": "26127749-a31f-40e6-b9cb-6792ce9fc8cb", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "shape: (1, 13)
termestimatestd_errorstatisticp_values_valueconf_lowconf_highstatistic_noninfstatistic_nonsupp_value_noninfp_value_nonsupp_value_equiv
strf64f64f64f64f64f64f64f64f64f64f64f64
"b2-b1=0"0.6372690.1880083.3895850.000710.4803860.268781.0057573.336395-1.9293410.0004240.0268440.026844
" ], "text/plain": [ "shape: (1, 8)\n", "┌─────────┬──────────┬───────────┬──────┬─────────┬──────┬───────┬───────┐\n", "│ Term ┆ Estimate ┆ Std.Error ┆ z ┆ P(>|z|) ┆ S ┆ 2.5% ┆ 97.5% │\n", "│ --- ┆ --- ┆ --- ┆ --- ┆ --- ┆ --- ┆ --- ┆ --- │\n", "│ str ┆ str ┆ str ┆ str ┆ str ┆ str ┆ str ┆ str │\n", "╞═════════╪══════════╪═══════════╪══════╪═════════╪══════╪═══════╪═══════╡\n", "│ b2-b1=0 ┆ 0.637 ┆ 0.188 ┆ 3.39 ┆ 0.0007 ┆ 10.5 ┆ 0.269 ┆ 1.01 │\n", "└─────────┴──────────┴───────────┴──────┴─────────┴──────┴───────┴───────┘\n", "\n", "Columns: term, estimate, std_error, statistic, p_value, s_value, conf_low, conf_high, statistic_noninf, statistic_nonsup, p_value_noninf, p_value_nonsup, p_value_equiv" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Equivalence hypothesis\n", "eq = me.predictions(looks_mod,\n", " newdata=dg,\n", " by='looks',\n", " hypothesis='b2 - b1 = 0',\n", " equivalence=[0.01, 1.])\n", "eq" ] }, { "cell_type": "markdown", "id": "f494e0f4-7f1b-4e66-a5a5-a9da38cc810b", "metadata": {}, "source": [ "This *is* significant, which indicates the effect *is within the interval*, and thus you can conclude the end result is the same as what you presently earn:" ] }, { "cell_type": "code", "execution_count": 16, "id": "52ed376c-cd4d-4776-af15-a7d0a9c6cd59", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlEAAAHHCAYAAACfqw0dAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA5nElEQVR4nO3de1jUZf7/8dcwjgKelRABD+UBMcnDmmSWeEJ/mZUaua6HNrfSb+rm2rey0vK7pFnpqil52E7mYXN1c2vJqzRRv2m6Jm4eQwg8hCdAARUHdID5/eGXKZLD8JlRBng+rstrdz73Pff9nq4LeXl/7rk/JrvdbhcAAAAqxKuyCwAAAKiKCFEAAAAGEKIAAAAMIEQBAAAYQIgCAAAwgBAFAABgACEKAADAgFqVXUB1dujQIdlsNnl5ealOnTqVXQ4AAHDC1atXVVhYKIvForCwsFL7EaJuIpvNJrvdroKCAlmt1souBwAAVIDNZiuznRB1E3l5eamgoEAmk0k+Pj6VXQ4AAHBCbm6u7Ha7vLzK3vVEiLqJ6tSpI6vVKh8fH4WGhlZ2OQAAwAkJCQmyWq3lbsVhYzkAAIABhCgAAAADCFEAAAAGEKIAAAAMIEQBAAAYQIgCAAAwgBAFAABgACEKAADAAEIUAACAAVUyRJ04cUJdunRRdHR0hd+blpammTNnKjIyUmFhYerbt69ef/11ZWZm3oRKAQBAdVXlQtT58+c1ceJE5ebmVvi9qampevTRR7V27Vp5e3urb9++MpvNWr16tYYNG6Zz587dhIoBAEB1VKVCVEJCgkaNGqWUlBRD73/ppZeUkZGhSZMmKTY2VosWLdKmTZs0cuRInTt3TjNnznRzxQAAoLqqEiHq4sWLmjt3rkaMGKGTJ08qODi4wmPEx8crPj5erVu31uTJkx3XzWazZsyYocDAQG3fvl3JycnuLB0AAFRTVSJErVy5Uu+//76aNGmipUuXaujQoRUeY+vWrZKk/v37y8ur+Me2WCzq16+fJCkuLs7legEAQPVXq7ILcEZAQICmTZumUaNGydvbW0eOHKnwGElJSZKk9u3bl9jetm1bSdLRo0eNF1qK/IJCZeXkuX1cAJXjqq1AH+064Xg97t7WqmMxV15BANwqv6DQqX5VIkQ99thjLo+Rnp4uSWrWrFmJ7f7+/sX6uVPCuct6aDUrXEB1tXzHicouAYAbvT2gqdo0tpTbr0rcznMHq9UqSfLx8Smx3dvbu1g/AACAstSYEGU2X19qN5lMZfaz2+23ohwAAFDFVYnbee5Qt25dSSr1fKm8vOt7lkpbqXJFx+YNdGDmb9w+LoDKcSnPpvvf2uZ4vWNaXzXwLn/pH0DVkHrsR13NK/88yhoTovz9/XXkyBFlZGSU2F60F6pob5Q7mb1MaujDX7BAddXA28LPOFCNnPEq+65VkRpzOy8kJESSSj0Hquh6UT8AAICy1JgQFRERIUn6+uuvVVhY/KuLNpvNcT5U3759b3ltAACg6ql2IcpmsyklJUUpKSmy2WyO6926dVNYWJhSUlI0f/58xwbygoICzZ49W2fPnlXv3r0VGhpaWaUDAIAqpNrtiUpLS9PgwYMlXT99/JePiJkzZ47GjBmj9957T3FxcWrXrp0SEhL0008/KSgoSLNmzaqssgEAQBVT7VaiytKuXTtt2LBBw4cP1+XLl7Vt2/Vv14wdO1br1q0r9SBOAACAXzPZORjppklISJDVapWvry+3CYFq5GKuTZ3/vNnx+sDMgXw7D6hGnP39XaNWogAAANyFEAUAAGAAIQoAAMAAQhQAAIABhCgAAAADCFEAAAAGEKIAAAAMIEQBAAAYQIgCAAAwgBAFAABgACEKAADAAEIUAACAAYQoAAAAAwhRAAAABhCiAAAADCBEAQAAGECIAgAAMIAQBQAAYAAhCgAAwABCFAAAgAGEKAAAAAMIUQAAAAYQogAAAAwgRAEAABhAiAIAADCAEAUAAGAAIQoAAMAAQhQAAIABhCgAAAADCFEAAAAGEKIAAAAMIEQBAAAYQIgCAAAwgBAFAABgACEKAADAAEIUAACAAYQoAAAAAwhRAAAABhCiAAAADCBEAQAAGECIAgAAMIAQBQAAYAAhCgAAwABCFAAAgAGEKAAAAAMIUQAAAAYQogAAAAwgRAEAABhAiAIAADCAEAUAAGBArcouwFnHjx/Xu+++q3379unChQsKCAjQAw88oAkTJsjX17dCY+3fv1/Lly/Xf/7zH125ckVNmzbVvffeq2eeeUYtW7a8SZ8AAABUJ1ViJergwYMaPny4YmNj5efnpz59+shqtWrZsmUaOXKkcnJynB7ryy+/1KhRo7R161YFBQWpT58+ql27tjZs2KChQ4fq4MGDN/GTAACA6sLjQ1R+fr6ee+45Wa1WzZ49W+vXr9eiRYu0ZcsW9evXT4mJiZo/f75TY+Xm5uq1115TYWGh/vKXv2jDhg2KiYnRV199pSeeeEJXrlzR9OnTb/InAgAA1YHHh6iNGzcqNTVVPXv2VFRUlOO6t7e33njjDfn6+mrdunW6ePFiuWPt379fly5dUrt27TRkyBDHdbPZrKlTp8psNispKUmZmZk35bMAAIDqw+ND1NatWyVJkZGRN7Q1btxY4eHhstls2rFjR7ljmc1mSVJmZqauXbtWrC0rK0sFBQWyWCyqV6+eGyoHAADVmceHqKSkJElSSEhIie1t27aVJB09erTcscLCwtSkSROdP39e//3f/61jx44pLy9PBw4c0KRJkyRJTzzxhGrXru2m6gEAQHXl8d/OS09PlyQ1a9asxHZ/f/9i/cri4+OjmJgYTZ06VZs3b9bmzZsdbd7e3oqOjtZvf/tbN1RdnN1uV0FBgdvHBVA5Cn/181xYUKCCAo//NykAJ9ntdqf6eXyIslqtkq6HnJIUXS/qV57bb79dDz/8sD788EN16NBBAQEBSkpKUmpqqlasWKGOHTsqLCzMPcX/n9zcXO3fv9+tYwKoPFeuFRZ7fejQIdWtTYgCahqPD1Fms1mFhYUymUxl9nMmNZ49e1ZjxoxRVlaWPvjgA/Xs2dPx3o8//lhz5szRuHHjFBsbq+bNm7ulfgAAUD15fIiqW7eusrOzlZubW2J7Xl6epOu36sqzYMECnTp1Si+//LIjQEmSyWTSE088ocOHDys2NlYrVqzQyy+/7J4P8H+1lbanC0DVcynXJn0e53gdFhamBj6WSqwIgDslJiaWmjt+yeNDlL+/v7Kzs5WRkaEWLVrc0F60F6pob1RZdu3aJUm6//77S2zv06ePYmNjdfjwYRcqvpHJZHJ8MxBA1edlLvzVazM/40A1Ut7dryIefxO/aAUnOTm5xPai686s9BSdJVWrVsnZsegvwfz8/ArXCQAAahaPD1ERERGSpE2bNt3QlpWVpT179shisahXr17ljtWmTRtJP5899Ws7d+6UJIWGhhotFwAA1BAeH6IiIyMVGBionTt3as2aNY7reXl5mj59uqxWq6KiouTn5+dos9lsSklJUUpKimw2m+P66NGjJUkxMTHau3dvsXn+8Y9/6NNPP5XFYtGYMWNu8qcCAABVncfvifL29tabb76p8ePHKzo6Wp9++qmCg4P1/fffKz09XR07dtTzzz9f7D1paWkaPHiwJCkuLk7BwcGSpMcee0yHDh3S3//+d40ZM0ZhYWEKCAhQcnKyjh8/LovFotmzZzsO8AQAACiNx4coSQoPD9f69esVExOj7777TsnJyQoODlZUVJSefPLJCj2mJTo6Wr1799Ynn3yiw4cPKyEhQY0bN9aQIUP01FNPcSsPAAA4pUqEKElq3769Fi1a5FTf4OBgJSYmlto+YMAADRgwwF2lAQCAGsjj90QBAAB4IkIUAACAAYQoAAAAAwhRAAAABhCiAAAADCBEAQAAGECIAgAAMIAQBQAAYAAhCgAAwABCFAAAgAGEKAAAAAMIUQAAAAYQogAAAAwgRAEAABhAiAIAADCAEAUAAGAAIQoAAMAAQhQAAIABhCgAAAADCFEAAAAGEKIAAAAMIEQBAAAYQIgCAAAwgBAFAABgACEKAADAAEIUAACAAYQoAAAAAwhRAAAABhCiAAAADCBEAQAAGECIAgAAMIAQBQAAYAAhCgAAwABCFAAAgAGEKAAAAAMIUQAAAAYQogAAAAwgRAFABaRfytOq3SeKXVu1+4TSL+VVTkEAKk2tyi4AAKqCQ6cuatk3Kdp0+JzyC+3F2uZtTtLCLT9qUKcA/VfvNgoLblhJVQK4lQhRAFCO1f8+qdc+P6xfZadi8gvt2njwrL48dFbRj3TSmHta3boCAVQKQhQAlGH1v09qxmeHne5faJdmfHZYJpM0OpwgBVRn7IkCgFIcOnVRr33ufID6pVc/O6xDpy66uSIAnoQQBQClWPZNSpm38MpSaJeWf5Pi3oIAeBRCFACUIP1SnjYdPufSGF8dPse39oBqjBAFACX414EzN3wLr6LyC+3614EzbqoIgKdx28byixcvavfu3Tp27JguX76sadOm6erVqzpw4IB69OjhrmkA4JZIychx0zhX3DIOAM/jcoiy2+1avHixPvroI+Xl/bxsPW3aNJ06dUq///3v1aVLF7377rtq0qSJq9MBwC1x5WqBm8bJd8s4ADyPy7fzXnzxRS1dulS5ublq1KiRfHx8HG3Z2dmy2+3av3+/xo4dq9zcXFenA4Bbom4ds5vG4SQZoLpyKURt3rxZsbGxatKkid577z3t3r1bHTp0cLT/5je/0Zo1a9S0aVMdO3ZMK1eudLlgALgV2txWz03j1HXLOAA8j0sh6u9//7tMJpP+8pe/6P777y+xz29+8xstXLhQdrtdmzZtcmU6ALhlHu4cqFpeJpfGqOVl0sOdA91UEQBP41KIOnz4sJo3b6577rmnzH7du3dXUFCQTpw44cp0AHDL+Dfw1qBOAS6N8f86Bci/gbebKgLgaVwKUVarVY0aNXKqb5MmTZSfb3yD5fHjx/X888+rb9++uuuuuzRw4EAtWLBAVqu1wmNZrVbFxMTooYceUufOndW1a1eNGjVKX331leH6AFQ//9W7jYwuRnmZpAm927i3IAAexaUQ5efnp5MnT8puL/ssFZvNphMnTsjPz8/QPAcPHtTw4cMVGxsrPz8/9enTR1arVcuWLdPIkSOVk+P8V5HPnz+vxx57TIsXL1ZmZqbuu+8+hYSE6D//+Y+mTJmiFStWGKoRQPUTFtxQ0Y90MvTe14d2UlhwQzdXBMCTuBSievToIavVqr/97W9l9vv44491+fJlde/evcJz5Ofn67nnnpPVatXs2bO1fv16LVq0SFu2bFG/fv2UmJio+fPnOz3e9OnTlZycrEGDBmnr1q169913tXbtWn344YeyWCx6++23deYMh+MBuG7MPa00e1gnp1ekvEzS7GGdePgwUAO4FKLGjRsnLy8vvfXWW1q5cqWysrKKtV+4cEELFy7U/Pnz5eXlpTFjxlR4jo0bNyo1NVU9e/ZUVFSU47q3t7feeOMN+fr6at26dbp4sfwHfR48eFDbt29Xy5YtNXfuXNWpU8fRdu+992r48OEKCAjQ4cPGHjgKoHoaHd5Kn0+6T0Pual7qZvNaXiYNuau5Pp90HwEKqCFcOsCkQ4cOeuWVVzRr1izNmTNHc+bMcbT17NlT2dnZkq4fyDllyhTdddddFZ5j69atkqTIyMgb2ho3bqzw8HBt27ZNO3bs0JAhQ8oc68svv5QkjR07tliAKhIdHV3h+gDUDGHBDRUzqpvSL+VpXXyq5m1OcrQ9P7C9RtzdQv712UQO1CQunwI3evRotWjRQvPnz9fRo0cd14tWpVq1aqUpU6Zo8ODBhsZPSrr+F1VISEiJ7W3bttW2bdt09OjRckNU0QpTly5dZLVatWnTJh06dEgFBQUKCwvTkCFD5O3t/r8E7Xa7Cgrcc/oxgMrVtK5Fo3u0KBaiRvdooQY+Fn7OgWqivL3eRdxylG7v3r3Vu3dvnT59Wj/++KMuX74sHx8f3X777WrTxrVvp6Snp0uSmjVrVmK7v79/sX5lKTpiISsrS0OGDNHp06cdbWvXrtXSpUu1bNkytWvXzqWafy03N1f79+9365gAKs+Va4XFXh86dEh1a/M8d6CmcSlE7d27V/Xr13ecUh4UFKSgoKAS++7cuVMnT57U6NGjKzRH0REGpa0QFV135qiDom/xPf/88woMDNTq1asVGhqqU6dOad68edqxY4eefvppffHFF6pXzz2nFQMAgOrJpRA1duxYde/eXatXry6378KFCw2FKLPZrMLCQplMZX81xpmlt6tXr0qSateurZUrV6phw+tfP+7QoYOWLVumYcOGKSkpSWvXrtVTTz1VoTrL4uPjU+rtSABVz6Vcm/R5nON1WFiYGvhYKrEiAO6UmJjo1PN+nQ5ROTk5N3z7TpLy8vKUmppa6vvsdrtOnz6tY8eOOX2P8Zfq1q2r7OzsUj9MXl6eJBV78HFpvL29deXKFQ0dOtQRoIrUqlVLI0eOVHR0tHbv3u3WEGUymWQ2u+dhpgAqn5e58FevzfyMA9VIeQs3RZwOUVeuXNFDDz3kWM0pmuTIkSMaOHCgU2N06dLF2ekc/P39lZ2drYyMDLVo0eKG9qK9UEV7o8ri5+enK1euKDg4uMT2ouslhUUAAIBfcnonZLNmzTRu3DjZ7XbHH0nFXpf1p3nz5poxY0aFCyy6DZacnFxie9F1Z26XFfVJS0srsT0jI0PS9UfUAAAAlKVCe6ImTpzoOPDSbrdrwIABCgsL08KFC0t9j5eXl3x9fW+4feasiIgIxcbGatOmTRoxYkSxtqysLO3Zs0cWi0W9evUqd6w+ffpo8+bN+vLLLzVp0iRZLMX3MHzzzTeSrp/EDgAAUJYKfSfXYrE4voEXHBysu+++W507d3ZcK+lP8+bNDQco6fohm4GBgdq5c6fWrFnjuJ6Xl6fp06fLarUqKiqq2HP5bDabUlJSlJKSIpvN5rg+ePBgBQcH68SJE4qOji7Wtn79em3atEkNGzbUo48+arheAABQM5jsRnZ732J79uzR+PHjlZeXpzvvvFPBwcH6/vvvlZ6ero4dO2rVqlXFjiQ4deqU+vfvL0mKi4srtgfq8OHDeuqpp5SVlSV/f3917txZJ0+eVFJSkurUqaOFCxeqX79+bqk7ISFBVqtVvr6+Cg0NdcuYACrfxVybOv95s+P1gZkD1ZBv5wHVhrO/v91+OlxhYWGxP9euXdOlS5eUlJSk5cuXGxozPDxc69ev16BBg3TmzBlt375d9evX18SJE28IUOXp1KmTYmNjNXbsWNWuXVvbt29XVlaWHnzwQa1bt85tAQoAAFRvLq9EffPNN1q8eLGSkpJ07dq1cvsnJCS4Ml2VwkoUUD2xEgVUb87+/nbpsM3Dhw9r4sSJKigoKPcMqFq1aqlr166uTAcAAOAxXApRK1euVH5+vtq1a6ennnpK3t7emjJligYNGqTf/va3OnfunDZs2KD4+Hh1795dK1ascFPZAAAAlculEBUfHy+z2axFixbp9ttvlyQ1b95cqampuvfeeyVJw4YN07PPPqstW7Zo48aNevDBB12vGgAAoJK5tLH8woULCgwMdAQo6fpz6H65P8pkMunll1+WJG3YsMGV6QAAADyGy9/Oa9SoUbHXrVu3VkFBgY4fP+64FhgYqFatWikxMdHV6QAAADyCSyGqadOmjkelFCk6k+nXj2kpepAwAABAdeBSiOrUqZPS0tK0Y8cOx7U77rhDdrtde/fudVzLy8vTTz/95NLJ5QAAAJ7EpRD1yCOPyG6369lnn9XcuXOVn5+vrl27qmHDhvrHP/6hzz77TElJSZoxY4YuX75cbO8UAABAVeZSiOrfv78efPBB5ebm6uOPP5bZbJa3t7d+97vfKT8/Xy+//LIeeeQRbdy4USaTSY8//ri76gYAAKhULh1xIEl/+ctf1KtXL3377bcymUySpMmTJystLU2fffaZ7Ha7zGazxo0bp4EDB7pcMAAAgCdwOURJ0vDhwzV8+PCfB61VS3PmzNHUqVN15swZtWzZUk2aNHHHVAAAAB7BLSGqNP7+/vL397+ZUwAAAFQKt4Soa9eu6dixY8rJySn3GXp33323O6YEAACoVC6HqAULFmjlypXKy8srt6/JZNIPP/zg6pQAAACVzqUQtWLFCi1fvtzp/uWtUgEAAFQVLoWodevWyWQyafjw4Zo0aZL8/f1Vq9ZN3WYFAADgEVxKPKmpqWratKlef/11eXm5/Bg+AACAKsOl5FOvXj3ddtttBCgAAFDjuJR+7r77bh0/flxXrlxxVz0AAABVgkshauLEiSooKNCf//xnNo0DAIAaxaU9UR06dNA777yjyZMna//+/erVq5eaNGniePxLSSZPnuzKlAAAAB7BpRCVk5Ojv/71ryosLFRqaqrWrl1b7nsIUQAAoDpwKUQtWLBA+/fvlyTddtttCgwMlMVicUddAAAAHs2lEBUXFyeTyaQ33nhDw4YNc1dNAAAAHs+ljeWZmZlq1aoVAQoAANQ4LoUof39/mc1md9UCAABQZbgUogYNGqRjx47pyJEj7qoHAACgSnApRD3zzDNq2bKlnnnmGX311VccugkAAGoMlzaWR0dHq0WLFtq5c6emTp0qk8mk+vXry8fHp8T+JpNJ27Ztc2VKAAAAj+BSiPrXv/7l+P92u112u10XL17UxYsXS+xf1iGcAAAAVYlLIWrOnDnuqgMAAKBKcSlEcbQBAACoqVzaWA4AAFBTOb0StXv3bklSt27dVKdOnWLXKqJnz54Vfg8AAICncTpEjRs3Tl5eXtq4caNuv/12x7WKbBY3mUz64YcfKl4lAACAh6nQnqjCwsIbrtntdqffX5G+AAAAnszpEHX06FGnrgEAANQEt3Rj+YULF27ldAAAADeNSyGqf//+mjp1qlN9R44cyZEIAACg2nApRJ0+fVrp6enl9issLFRGRoaysrJcmQ4AAMBjOL0nKjk5WTNnzrzhelJSkkaPHl3q++x2u9LS0nTmzBkFBgYaqxIAAMDDOB2i2rZtK29vb3377beOayaTSZcvX9a+ffucGmPMmDEVrxAAAMADVeiIg1dffVVffPGF43VMTIwCAwM1fPjwUt9jMplUt25dhYaGKjw83HilAAAAHqRCIap169aaPHmy43VMTIyaN29e7BoAAEBN4NIDiOPi4hyPgAEAAKhJXApRQUFBpbbl5eVp165dKiwsVPfu3dWoUSNXpgIAAPAoLoUoSUpLS9PSpUsVGBio8ePHS5JSUlI0btw4ZWRkSJJ8fHw0a9YsDR482NXpAAAAPIJLISozM1MjRoxQenq6+vTp47j+2muvKT093bGpPCcnRy+++KJCQkLUpk0bV2sGAACodC4dtvnxxx8rLS1NLVu21G9/+1tJ0smTJ7Vv3z6ZzWZ98sknio+P1/jx45Wfn68VK1a4o2YAAIBK51KI+uabb1SrVi198MEHjpWo7du3S5K6deumLl26SJL++Mc/qkGDBvr3v//tynQAAAAew6UQlZqaqtatWys4ONhxbdeuXTKZTLr33nsd1ywWi4KDg516RExpjh8/rueff159+/bVXXfdpYEDB2rBggWyWq2ufARJ0ltvvaWQkBAtXrzY5bEAAEDN4FKIKigoUO3atR2v8/PztXfvXklSjx49ivXNzc2VyWQyNM/Bgwc1fPhwxcbGys/PT3369JHVatWyZcs0cuRI5eTkGP4M3377rT766CPD7wcAADWTSyEqKChIp0+fls1mkyTt3btXVqtVdevWddzKk65/gy81NVXNmzev8Bz5+fl67rnnZLVaNXv2bK1fv16LFi3Sli1b1K9fPyUmJmr+/PmG6s/MzNS0adNkt9sNvR8AANRcLoWosLAwXbp0SfPmzdPRo0e1cOFCmUwmRUREyGw2S5IuXLigF154QQUFBerZs2eF59i4caNSU1PVs2dPRUVFOa57e3vrjTfekK+vr9atW6eLFy9WeOxXXnlFWVlZ6tatW4XfCwAAajaXQtTTTz8tb29vrVy5UsOGDdOBAwdkNpv19NNPS5Li4+MVERGhvXv3qn79+vrDH/5Q4Tm2bt0qSYqMjLyhrXHjxgoPD5fNZtOOHTsqNO6aNWu0bds2TZo0SZ06dapwXQAAoGZzKUTdcccd+vDDDxUWFqbatWurffv2Wrp0qTp06CBJ8vf3V35+vtq1a6dPPvmk2AZ0ZyUlJUmSQkJCSmxv27atJOno0aNOj/njjz/qrbfeUrdu3TRhwoQK1wQAAODyieVdu3bVunXrSmwLDg7WZ5995ghVRhR9o69Zs2Yltvv7+xfrV56rV6/queeek8Vi0dy5cx23HW8mu92ugoKCmz4PgFuj8Fc/z4UFBSoocOnfpAA8iLN7pV0OUUV27Nihbdu26dixY7p8+bI+/fRT5eTkaMuWLWrWrJkaN25saNyiIwy8vb1LbC+67uxRB2+//baSkpL01ltvGVoZMyI3N1f79++/JXMBuPmuXCss9vrQoUOqW5sQBdQ0LoeoCxcu6E9/+pPi4+MlXU9vRUcZnDlzRjExMVq1apX++te/qnPnzhUe32w2q7CwsNzjEZxJjdu3b9fq1as1ePBgDR06tMK1AAAAFHEpRF27dk1PPvmkjh49qnr16unee+/VgQMHHLfWvLy81KhRI2VnZ2vcuHGKjY1VUFBQheaoW7eusrOzlZubW2J7Xl6epOsPOS7L+fPn9fLLL6t58+b685//XKEaXOXj41Pqni4AVc+lXJv0eZzjdVhYmBr4WCqxIgDulJiYWGru+CWXQtSaNWt09OhRdenSRUuWLFGTJk00atQoR4hq3769tmzZoqefflr79+/XRx99pBkzZlRoDn9/f2VnZysjI0MtWrS4ob1orqK9UaVZsmSJMjMzFRoaqujo6GJtR44ckSRt3rxZJ0+eVJs2bfTMM89UqM6ymEymW7L3CsCt4WUu/NVrMz/jQDXi7OHgLoWojRs3ysvLS3PnzlWTJk1K7FOvXj3NmzdPgwYNqvAxBNL1b+UlJSUpOTm5xPOckpOTHf3KUrRnKiEhQQkJCSX2SUpKUlJSknr06OHWEAUAAKofl3ZCHjt2TG3atClxheiXgoKC1Lp1a509e7bCc0REREiSNm3adENbVlaW9uzZI4vFol69epU5zptvvqnExMQS/zz++OOSpMmTJysxMVGrVq2qcJ0AAKBmcSlEFRYWlt/p/1gsFkPL3ZGRkQoMDNTOnTu1Zs0ax/W8vDxNnz5dVqtVUVFR8vPzc7TZbDalpKQoJSXF8UgaAAAAd3Lpdl5QUJBOnDihnJwc1atXr9R+WVlZ+vHHH9W6desKz+Ht7a0333xT48ePV3R0tD799FMFBwfr+++/V3p6ujp27Kjnn3++2HvS0tI0ePBgSVJcXNwtO8oAAADUHC6tREVERMhms2nu3Lll9ps1a5YKCgp0//33G5onPDxc69ev16BBg3TmzBlt375d9evX18SJE7Vq1aoyAxwAAMDNYLI7eyxnCS5cuKDBgwfr0qVL6t+/vx566CHFxMQoOTlZX3zxhZKSkrRmzRrt27dPdevW1caNG0s9ebw6SkhIkNVqla+vr0JDQyu7HABucjHXps5/3ux4fWDmQDXkiAOg2nD297dLt/OaNm2qJUuWaOLEidqyZYvi4n4+N2XIkCGSrh+C6evrq/nz59eoAAUAAKo3l59T8Jvf/Eb/+te/9Pjjj6t58+ay2+2OP02bNlVUVJQ+++wz9e7d2x31AgAAeAS3PDuvWbNmeuWVV/TKK6/IarXq8uXL8vX1Vf369d0xPAAAgMdx2wOIi/j6+srX19fdwwIAAHgUHjsOAABgACEKAADAAEIUAACAAYQoAAAAAwhRAAAABhCiAAAADCBEAQAAGECIAgAAMIAQBQAAYAAhCgAAwABCFAAAgAGEKAAAAAMIUQAAAAYQogAAAAwgRAEAABhAiAIAADCAEAUAAGAAIQoAAMAAQhQAAIABhCgAAAADCFEAAAAGEKIAAAAMIEQBAAAYQIgCAAAwgBAFAABgACEKAADAAEIUAACAAYQoAAAAAwhRAAAABhCiAAAADCBEAQAAGFCrsgsAgKqmTi0vTenfrthrADUPIQoAKsjbYtbUyPaVXQaASsY/nwAAAAwgRAEAABhAiAIAADCAEAUAAGAAIQoAAMAAQhQAAIABhCgAAAADCFEAAAAGEKIAAAAMIEQBAAAYQIgCAAAwgBAFAABgACEKAADAAEIUAACAAbUquwBnHT9+XO+++6727dunCxcuKCAgQA888IAmTJggX1/fCo21fft2rV69WocPH9bly5fVqFEjdevWTU899ZQ6d+58kz4BAACoTqrEStTBgwc1fPhwxcbGys/PT3369JHVatWyZcs0cuRI5eTkOD3W/PnzNWHCBO3cuVNBQUHq06ePGjRooM2bN+t3v/ud/vnPf97ETwIAAKoLj1+Jys/P13PPPSer1arZs2crKipKkpSXl6epU6dq69atmj9/vl577bVyx4qPj9fy5cvl4+Oj5cuXKzw83NG2du1azZw5U6+99prCw8MVGBh40z4TAACo+jx+JWrjxo1KTU1Vz549HQFKkry9vfXGG2/I19dX69at08WLF8sd6x//+Ick6amnnioWoCRp5MiRioiI0LVr17Rp0yb3fggAAFDteHyI2rp1qyQpMjLyhrbGjRsrPDxcNptNO3bsKHcsb29vtW/fXvfcc0+J7XfccYckKS0tzYWKAQBATeDxt/OSkpIkSSEhISW2t23bVtu2bdPRo0c1ZMiQMsf6n//5nzLbDxw4IEkKCAioeKFlsNvtKigocOuYAADg5rDb7U718/gQlZ6eLklq1qxZie3+/v7F+hm1detW/ec//5HFYilx1csVubm52r9/v1vHBAAAlcvjb+dZrVZJ12/FlaToelE/IxITE/Xyyy9Lur5fKigoyPBYAACgZvD4lSiz2azCwkKZTKYy+zm79PZrBw8e1Pjx45Wdna2+ffvq2WefNTROWXx8fEq9HQkAADxLYmKicnNzy+3n8SGqbt26ys7OLvXD5OXlSboeVCrqq6++0ksvvaTc3FwNGDBACxYskJeX+xfnTCaTzGaz28cFAADuV97CTRGPv51XtOcpIyOjxPaivVBF/Zz17rvv6k9/+pNyc3M1atQoLVq0SLVr13atWAAAUGN4fIgqug2WnJxcYnvRdWdvlxUWFuqll17SokWL5OXlpenTp2vmzJmsFAEAgArx+BAVEREhSSUegJmVlaU9e/bIYrGoV69eTo03Y8YM/fOf/5Svr6+WLl2qxx9/3K31AgCAmsHjQ1RkZKQCAwO1c+dOrVmzxnE9Ly9P06dPl9VqVVRUlPz8/BxtNptNKSkpSklJkc1mc1z/7LPP9Omnn6pWrVpasmSJI6ABAABUlMlu9Gttt9CePXs0fvx45eXl6c4771RwcLC+//57paenq2PHjlq1apXq1avn6H/q1Cn1799fkhQXF6fg4GAVFBSof//+Onv2rJo1a6YePXqUOt99992noUOHulx3QkKCrFarfH19FRoa6vJ4AADg5nP297fHfztPksLDw7V+/XrFxMTou+++U3JysoKDgxUVFaUnn3yyWIAqTWJios6ePSvp+mNdYmNjS+3boEEDt4QoAABQfVWJlaiqipUoAACqHmd/f3v8nigAAABPRIgCAAAwgBAFAABgACEKAADAAEIUAACAAYQoAAAAAwhRAAAABhCiAAAADCBEAQAAGECIAgAAMIAQBQAAYAAhCgAAwABCFAAAgAGEKAAAAAMIUQAAAAYQogAAAAwgRAEAABhAiAIAADCAEAUAAGAAIQoAAMAAQhQAAIABhCgAAAADCFEAAAAGEKIAAAAMIEQBAAAYQIgCAAAwgBAFAABgACEKAADAAEIUAACAAYQoAAAAAwhRAAAABhCiAAAADCBEAQAAGECIAgAAMIAQBQAAYAAhCgAAwABCFAAAgAGEKAAAAAMIUQAAAAYQogAAAAwgRAEAABhAiAIAADCAEAUAAGAAIQoAAMAAQhQAAIABhCgAAAADCFEAAAAGEKIAAAAMIEQBAAAYQIgCAAAwoFZlF+Cs48eP691339W+fft04cIFBQQE6IEHHtCECRPk6+tbobHS0tK0ZMkS7dq1S+fOnZOfn5/69eunSZMmqUmTJjfpEwAAgOqkSqxEHTx4UMOHD1dsbKz8/PzUp08fWa1WLVu2TCNHjlROTo7TY6WmpurRRx/V2rVr5e3trb59+8psNmv16tUaNmyYzp07dxM/CQAAqC48fiUqPz9fzz33nKxWq2bPnq2oqChJUl5enqZOnaqtW7dq/vz5eu2115wa76WXXlJGRoYmTZqkZ599VpJUUFCg6OhorV27VjNnztTy5cvd+hkKCu26mGtz65gAAODmKCi0O9XPZLfbnetZST7//HO9+OKL6tmzp1asWFGsLSsrS/369ZPNZtO3336rhg0bljlWfHy8Ro8erdatW+vLL7+Ul9fPC3E2m00DBw7UmTNntHHjRrVt29bl2hMSEmS1WpWSZdOLWy64PB4AALj53h7QVG0aW+Tr66vQ0NBS+3n87bytW7dKkiIjI29oa9y4scLDw2Wz2bRjxw6nx+rfv3+xACVJFotF/fr1kyTFxcW5WjYAAKjmPD5EJSUlSZJCQkJKbC9aMTp69KjTY7Vv397lsQAAQM3m8Xui0tPTJUnNmjUrsd3f379Yv1s1VkWEBtTX9zO6uHVMAABwc5w6kaJrV/PK7efxIcpqtUqSvL29S2wvul7Uz5mxfHx8XB6rIq5dzdOxxCNuHRMAAFQuj7+dZzabJUkmk6nMfs7sj3fnWAAAoGbz+JWounXrKjs7W7m5uSW25+VdX24rbXXp12NJcstYFeHj41Pqni4AAOBZEhMTS80Kv+TxIcrf31/Z2dnKyMhQixYtbmgv2r9UtJ+pvLGOHDmijIyMEtsrMlZFmEwmxyoYAADwbOXdsSri8bfzilZwkpOTS2wvuu7MSo87xwIAADWbx4eoiIgISdKmTZtuaMvKytKePXtksVjUq1cvp8f6+uuvVVhYWKzNZrM5zofq27evq2UDAIBqzuNDVGRkpAIDA7Vz506tWbPGcT0vL0/Tp0+X1WpVVFSU/Pz8HG02m00pKSlKSUmRzfbz41a6deumsLAwpaSkaP78+Y4N5AUFBZo9e7bOnj2r3r17l3k6KQAAgFQFHvsiSXv27NH48eOVl5enO++8U8HBwfr++++Vnp6ujh07atWqVapXr56j/6lTp9S/f39J108fDw4OdrT9+OOPGjNmjLKzs3XHHXeoXbt2SkhI0E8//aSgoCB98sknpZ4jVVFFj30p79h4AADgOZz9/e3xK1GSFB4ervXr12vQoEE6c+aMtm/frvr162vixIk3BKjytGvXThs2bNDw4cN1+fJlbdu2TZI0duxYrVu3zm0BCgAAVG9VYiWqqmIlCgCAqqdarUQBAAB4GkIUAACAAYQoAAAAAzz+xPKq7OrVq5KuP2YmISGhkqsBAADOKHrkS9Hv8dIQom6iogM97Xa7rFZrJVcDAAAq4tcHc/8aIeomslgsstls8vLyUp06dSq7HAAA4ISrV6+qsLBQFoulzH4ccQAAAGAAG8sBAAAMIEQBAAAYQIgCAAAwgBAFAABgACEKAADAAEIUAACAAYQoAAAAAwhRAAAABhCiAAAADCBEAQAAGECIAgAAMIAQBQAAYAAhCgAAwABCFAAAgAGEKAAAAAMIUQA80uLFixUSEqLo6Gin35Obm6slS5bokUceUZcuXXTXXXfpgQce0Lx585SdnV2h+ffs2aOQkBANGTKkgpUXd+LECXXp0qVCn8NZaWlpmjlzpiIjIxUWFqa+ffvq9ddfV2ZmptvnAnAjQhSAaiE7O1sjRozQO++8o9OnT6tr167q0aOHMjMz9d5772n48OE6d+7cLa3p/PnzmjhxonJzc90+dmpqqh599FGtXbtW3t7e6tu3r8xms1avXq1hw4bd8s8K1ESEKADVwty5c5WUlKQePXpo8+bN+uijj/T+++/r66+/1v3336/Tp0/r1VdfvWX1JCQkaNSoUUpJSbkp47/00kvKyMjQpEmTFBsbq0WLFmnTpk0aOXKkzp07p5kzZ96UeQH8jBAFoMrLy8vTF198IUl688031aRJE0dbgwYN9NZbb8lkMmnHjh3Kysq6qbVcvHhRc+fO1YgRI3Ty5EkFBwe7fY74+HjFx8erdevWmjx5suO62WzWjBkzFBgYqO3btys5OdntcwP4Wa3KLgAAyrNr1y4tXrxYP/zwg7y9vXX33XdrwoQJCgsLkyRduHBBd955p+x2u4KCgm54f9OmTdWwYUNlZ2crPT1djRs3rtD8P/30k+bNm6fdu3crPz9foaGhGjNmjAYPHnxD35UrV+r9999XQECAZs6cqSNHjigmJqbUsUNCQpyqISgoSFu3bpUkx//2799fXl7F/y1ssVjUr18/rV69WnFxcWrbtq2zHxNABRGiAHi0Xbt26ZNPPlFAQIAiIiKUmpqqr7/+Wtu2bdPChQsVGRmpoKAg/e1vfyt1jJMnTyo7O1teXl5q1qxZhea/cOGCRowYIZvNpvDwcF25ckV79+7Vvn37dOjQIU2bNq1Y/4CAAE2bNk2jRo2St7e3jhw5Uub4Dz30kFN1/HJ1LSkpSZLUvn37EvsWBaejR486NTYAYwhRADza8ePHNWzYMEVHR6t27dqSpDVr1ig6OlqvvPKKunfvXu7K0rx58yRJ9957rxo1alSh+TMzMxUaGqoPPvhATZs2lSTt379fTz75pD788EP17t1bPXv2dPR/7LHHKjR+UW0VkZ6eLkmlBkJ/f/9i/QDcHOyJAuDRmjZtqldffdURoCRp9OjRioiI0KVLl/T555+X+f5ly5Zp8+bN8vb21osvvmiohujoaEeAkqQuXbromWeekXT99t2tZrVaJUk+Pj4ltnt7exfrB+DmIEQB8GiRkZGqW7fuDdcHDBggSfr3v/9d6nsXLVqkBQsWyMvLS2+88YbT+49+6Y477tBdd91V6vzfffddhcd0ldlsliSZTKYy+9nt9ltRDlBjcTsPgEdr0aJFidebN28uqeRbVteuXdOMGTP0+eefq1atWpozZ44efPDBYn2WLl1a4vEDAwcO1MCBA8udPzAwUJKUk5OjnJwc1atXz7kP9CtGNpYXhcrSzp/Ky8uTVPpKFQD3IEQB8Gh16tQp8XrRKovFYil2/cKFC5o0aZK+//571atXT++8847uu+++G96/a9euEleRWrVqVSxElTb/L/3yVmNFGdlY7u/vryNHjigjI6PEvkXBsmhvFICbgxAFwKOlpaWVeP3UqVOSfl6Rkq4fRfDEE0/o9OnTCg4O1rJly9SuXbsS379q1SqX5k9NTZV0fc+WKyHKyMbykJAQbdu2rdRzoIquG7l9CcB57IkC4NF27NhR4t6eL7/8UpJ0zz33SLoedh5//HGdPn1ad911l9atW1dqgKqIH374QefPny93/lspIiJCkvT111+rsLCwWJvNZlNcXJwkqW/fvre8NqAmIUQB8GhHjx7VwoULiwWp5cuX67vvvlOzZs308MMPS5JeeOEFnT17ViEhIVqxYkWxb9O5wmaz6cUXX9SVK1cc13bv3q333ntPZrNZf/jDH9wyT0V069ZNYWFhSklJ0fz58x3/bQoKCjR79mydPXtWvXv3Vmho6C2vDahJuJ0HwKN17dpVy5Yt06ZNmxQSEqLk5GQlJyerfv36WrRokXx9ffXtt99qz549kqT69euX+dy4KVOmlLpZvCTt27fX/v37NWDAAHXv3l1ZWVmKj4+XJL366qvq1KmTax/QoDlz5mjMmDF67733FBcXp3bt2ikhIUE//fSTgoKCNGvWrEqpC6hJCFEAPNqQIUM0btw4LV++XFu3blW9evX0yCOP6I9//KMjDG3bts3RvyjglOb3v/99hULU7bffrrlz52ru3LnauXOnTCaTevbsqQkTJlTKrbwi7dq104YNGxQTE6MdO3Zo27ZtCggI0NixY/Vf//Vf8vPzq7TagJrCZOcgEQAAgApjTxQAAIABhCgAAAADCFEAAAAGEKIAAAAMIEQBAAAYQIgCAAAwgBAFAABgACEKAADAAEIUAACAAYQoANVebm6uTp06VdllAKhmCFEAqrXY2FgNGjRIu3fvruxSAFQzhCgA1dqCBQuUlpZW2WUAqIYIUQAAAAYQogAAAAwgRAGolhYvXqyQkBCdPn1akjRjxgyFhIRo8eLFjj7nz5/X22+/rcGDB6tz587q2rWrHn30UX344Ye6evVqqWPOmzdPW7Zs0aBBg9SpUyf169dPGzdu1KlTpxQSEqLevXursLBQa9as0dChQ9W5c2fdc889mjRpklJSUiRJmZmZev3119WnTx916tRJERERio6O1uXLl2/NfyAALqtV2QUAwM3QvHlzdevWTYcPH9a1a9fUqlUrNW3aVM2bN5ck7du3TxMnTlR2drYsFotat24tu92uI0eO6PDhw/r888/1/vvv67bbbrth7L179+rDDz9Uw4YN1aZNG6WkpCg0NNTRXlhYqClTpmjz5s1q1qyZWrVqpWPHjmnLli3au3evli9frilTpigjI0OtWrVSYGCgTp48qTVr1uiHH37QJ598IpPJdMv+WwEwxmS32+2VXQQA3Cz9+vXT6dOnNWvWLD322GOSpLS0ND388MPKzs7WiBEj9MILL6hBgwaSpJ9++knPP/+8Dhw4oO7du2vNmjWOsRYvXqyYmBhJUmRkpObPn6/atWsrMzNTTZo00alTp9S/f39JUq1atTRr1iwNHTpUJpNJSUlJGjFihHJzc+Xl5aUOHTpowYIFat26tSTp008/1SuvvCJJWrVqlXr06HGr/hMBMIjbeQBqnA8++EDZ2dnq16+fXn/9dUeAkqSWLVtqyZIlqlevnuLj4/W///u/JY4xbdo01a5dW5LUpEmTG9qjoqI0bNgwx4pS+/btHQHLbrfrnXfecQQoSXr00UcVFBQkSfrhhx/c8jkB3FyEKAA1zpYtWyRJDz/8cIntfn5+6tWrlyRp27ZtN7TfdtttatGiRZlz9OnT54ZrRSHp9ttvV8uWLW9o9/f3lyTl5OSUOTYAz8CeKAA1ypUrVxybzZcsWaKVK1eW2K+oz7Fjx25oKwo7ZSnae/VLFotFUskrV79sZ5cFUDUQogDUKL9c5UlKSiq3f0nflqtTp0657/Px8Sm1zcuLmwBAdUCIAlCj/DLcxMbGqn379pVYDYCqjH8OAahRGjRoID8/P0lScnJyqf0SExOVkJCgixcv3qrSAFQxhCgA1VrRt+N+uc+oaNP36tWrVVhYeMN7Ll++rN///vcaOnSoPv7441tSJ4CqhxAFoFrz9fWV9PNGcUkaP368fH19tW/fPr3wwgvKzMx0tJ0+fVrjx49XVlaW6tevr9GjR9/ymgFUDeyJAlCtdezYUUlJSXr//ff1zTffKDIyUhMnTtTChQs1depUffHFF9q0aZPatm0rm82mEydOKD8/X76+vvrrX/+qpk2bVvZHAOChCFEAqrVp06YpNzdXu3bt0rFjxxzProuIiNDGjRu1YsUK7dixQ8ePH1dBQYGCgoLUq1cv/eEPfyj3LCgANRuPfQEAADCAPVEAAAAGEKIAAAAMIEQBAAAYQIgCAAAwgBAFAABgACEKAADAAEIUAACAAYQoAAAAAwhRAAAABhCiAAAADCBEAQAAGECIAgAAMIAQBQAAYMD/B3hKnyFjtN3rAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Plot\n", "ax = sns.pointplot(data=eq, x='term', y='estimate')\n", "ax.vlines(eq['term'], eq['conf_low'], eq['conf_high'])\n", "ax.axhline(0.01)\n", "ax.axhline(1)" ] }, { "cell_type": "markdown", "id": "c6643063-27e2-40d3-9633-9189b58e528e", "metadata": {}, "source": [ "To reiterate, this test is significant, which means you can say the effect is *unintersting, null, or too small to be of interest*. If the test was *not* significant, you would conclude that the intervention on your looks could possibly yield effects outside of this interval and thus *could be of interest*." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.2" } }, "nbformat": 4, "nbformat_minor": 5 }